Re: [R] how to solve the step halving factor problems in gnls and nls

2005-07-19 Thread Dieter Menne
Yimeng Lu yl2058 at columbia.edu writes:

 Could you give me some advice in 
 solving the problem of such error message from gnls and nls?
 ## begin error message
 
 Problem in gnls(y1 ~ glogit4(b, c, m, t, x), data.frame(x..: Step halving 
 factor reduced below minimum in NLS step 
 

Try to set nlsTol in the optional control parameter (gnlsControl) to a larger 
value, e.g. 0.1 instead of default 0.001, but be sure to check your results 
make sense. If this should help, you can try intermediate values.

Dieter Menne

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


Re: [R] how to get dissimilarity matrix

2005-07-19 Thread Martin Maechler
 Baoqiang == Baoqiang Cao [EMAIL PROTECTED]
 on Mon, 18 Jul 2005 15:02:05 -0400 writes:

Baoqiang Hello All, I'm learning R. Just wonder, any
Baoqiang package or function that I can use to get the
Baoqiang dissimilarity matrix? Thanks.

Yes,
learn to use help.search()  {also read the docu :  ?help.search}

help.search(dissimilarity)

and find daisy() in recommended package 'cluster'.
There's also  dist() in 'stats' which is a bit less versatile.

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


Re: [R] New functions supporting GIF file format in R

2005-07-19 Thread Martin Maechler
 JarekT == Tuszynski, Jaroslaw W [EMAIL PROTECTED]
 on Mon, 18 Jul 2005 16:00:43 -0400 writes:

JarekT Hi, A minor announcement. I just added two functions
JarekT for reading and writing GIF files to my caTools
JarekT package. Input and output is in the form of standard
JarekT R matrices or arrays, and standard R color-maps
JarekT (palettes). The functions can read and write both
JarekT regular GIF images, as well as, multi-frame animated
JarekT GIFs. Most of the work is done in C level code
JarekT (included), so functions do not use any external
JarekT libraries.

 

JarekT For more info and examples go to
JarekT http://cran.r-project.org/doc/packages/caTools.pdf
JarekT http://cran.r-project.org/doc/packages/caTools.pdf
JarekT and click GIF.

Wouldn't it make sense to donate these to the 'pixmap' package
which is dedicated to such objects and has been in place for a
very long time?

Regards,
Martin

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


Re: [R] listing datasets from all my packages

2005-07-19 Thread Prof Brian Ripley
I did manage eventually to reproduce this via

data(package=makecdfenv)

which is a BioC package with a non-standard use of the data directory.
(There are no objects in the single R file in the data directory, which 
from its name I suggest is a dummy and the data directory should be 
removed entirely.  The indexing code fails to handle this correctly.)

I had no difficulty on a system that has almost all CRAN packages 
installed.

On Mon, 18 Jul 2005, Elizabeth Purdom wrote:

 Hi,
 I am using R 2.1.0 on Windows XP and when I type data() to list the
 datasets in R, there is a helpful hint to type 'data(package =
 .packages(all.available = TRUE))' to see the datasets in all of the
 packages -- not just the active ones.

 However, when I do this, I get the following message:
  data(package = .packages(all.available = TRUE))
 Error in rbind(...) : number of columns of matrices must match (see arg 2)
 In addition: Warning messages:
 1: datasets have been moved from package 'base' to package 'datasets' in:
 data(package = .packages(all.available = TRUE))
 2: datasets have been moved from package 'stats' to package 'datasets' in:
 data(package = .packages(all.available = TRUE))

 I possibly have old libraries in my R libraries because I copy them forward
 and update them with new versions of R, rather than redownload them. Is
 there a way to fix this or do the same another way? (I saw something in
 archives about a problem similar to this with .packages(), but I got the
 impression it was fixed for 2.1.0)

 Thanks,
 Elizabeth

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


-- 
Brian D. Ripley,  [EMAIL PROTECTED]
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] definition of index.array and boot.return in the code for boot

2005-07-19 Thread Prof Brian Ripley
Or read the *sources* for package boot, and find the original code with 
all the comments.

On Mon, 18 Jul 2005, Spencer Graves wrote:

 Excellent question.  Try 'getAnywhere(index.array)'.  It's hidden
 in namespace:boot.  Ditto for boot.return.

 spencer graves

 Obrien, Josh wrote:

 Dear R friends,

 I am reading the code for the function boot in package:boot in an attempt to 
 learn how and where it implements the random resampling used by the 
 non-parametric bootstraps.

 The code contains two (apparent) functions - 'index.array'  and  
 'boot.return' - for which I can find no documentation, and which don't even 
 seem to exist anywhere on the search path.  What are they?

 Also, if the meanings of those two don't answer my larger question, could 
 you point me to the code that implements the random resampling?

 Thanks very much for your help,

 Josh O'Brien
 UC Davis

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

 -- 
 Spencer Graves, PhD
 Senior Development Engineer
 PDF Solutions, Inc.
 333 West San Carlos Street Suite 700
 San Jose, CA 95110, USA

 [EMAIL PROTECTED]
 www.pdf.com http://www.pdf.com
 Tel:  408-938-4420
 Fax: 408-280-7915

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


-- 
Brian D. Ripley,  [EMAIL PROTECTED]
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] help: how to change the column name of data.frame

2005-07-19 Thread Prof Brian Ripley
Just change the names!  E.g.

names(DF)[c(4,6)] - names(DF)[c(6,4)]

Strictly a data frame has names, not column names, hence the use the 
names() and names-() functions here.

On Tue, 19 Jul 2005, wu sz wrote:

 I have a data frame with 15 variables, and want to exchange the data
 of 4th column and 6th column. First I append a column in the data
 frame, copy the 4th column data there, then copy the 6th column data
 to 4th column, and copy the appended column data to 6th column, but
 the names of the 4th and 6th column are still unchanged. How can I
 exchange them?

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] how to change bar colours in plot.stl

2005-07-19 Thread Michael Townsley
My thanks to Achim Zeileis and Prof Ripley for their responses.  In a very 
short time I not only had an answer and solved my problem, but also learned 
something about R that I can employ in other situations.

Much appreciated,

MT

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


Re: [R] colnames

2005-07-19 Thread Gilbert Wu
Hi Adai,

Many Thanks for the examples.

I work for a financial institution. We are exploring R as a tool to implement 
our portfolio optimization strategies. Hence, R is still a new language to us.

The script I wrote tried to make a returns matrix from the daily return indices 
extracted from a SQL database. Please find below the output that produces the 
'X' prefix in the colnames. The reason to preserve the column names is that 
they are stock identifiers which are to be used by other sub systems rather 
than R.

I would welcome any suggestion to improve the script.


Regards,

Gilbert

 p.RIs2Returns -
+ function (RIm)
+ {
+ x-RIm[1:(nrow(RIm)-1), 1:ncol(RIm)]
+ y-RIm[2:nrow(RIm), 1:ncol(RIm)]
+ RReturns - (y/x -1)
+ RReturns
+ }
 
 
 channel-odbcConnect(ourSQLDB)
 result-sqlQuery(channel,paste(select * from equityRIs;))
 odbcClose(channel)
 result
   stockidsdate  dbPrice
1   899188 20050713  7.59500
2   899188 20050714  7.60500
3   899188 20050715  7.48000
4   899188 20050718  7.41500
5   902232 20050713 10.97000
6   902232 20050714 10.94000
7   902232 20050715 10.99000
8   902232 20050718 11.05000
9   901714 20050713 17.96999
10  901714 20050714 18.00999
11  901714 20050715 17.64999
12  901714 20050718 17.64000
13  28176U 20050713  5.19250
14  28176U 20050714  5.25000
15  28176U 20050715  5.25000
16  28176U 20050718  5.22500
17  15322M 20050713 11.44000
18  15322M 20050714 11.5
19  15322M 20050715 11.33000
20  15322M 20050718 11.27000
 r1-reshape(result, timevar=stockid, idvar=sdate, direction=wide)
 r1
 sdate dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
dbPrice.15322M
1 20050713  7.595  10.97   17.96999 5.1925  
11.44
2 20050714  7.605  10.94   18.00999 5.2500  
11.50
3 20050715  7.480  10.99   17.64999 5.2500  
11.33
4 20050718  7.415  11.05   17.64000 5.2250  
11.27
 #Set sdate as the rownames
 rownames(r1) -as.character(r1[1:nrow(r1),1:1])
 #Get rid of the first column
 r1 - r1[1:nrow(r1),2:ncol(r1)]
 r1
 dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
dbPrice.15322M
20050713  7.595  10.97   17.96999 5.1925  
11.44
20050714  7.605  10.94   18.00999 5.2500  
11.50
20050715  7.480  10.99   17.64999 5.2500  
11.33
20050718  7.415  11.05   17.64000 5.2250  
11.27
 colnames(r1) - as.character(sub([[:alnum:]]*\\.,, colnames(r1)))
 r1
 899188 902232   901714 28176U 15322M
20050713  7.595  10.97 17.96999 5.1925  11.44
20050714  7.605  10.94 18.00999 5.2500  11.50
20050715  7.480  10.99 17.64999 5.2500  11.33
20050718  7.415  11.05 17.64000 5.2250  11.27
 RRs-p.RIs2Returns(r1)
 RRs
  X899188  X902232  X901714  X28176U  X15322M
20050714  0.001316656 -0.002734731  0.002225933  0.011073664  0.005244755
20050715 -0.016436555  0.004570384 -0.019988906  0.0 -0.014782609
20050718 -0.008689840  0.005459509 -0.000566006 -0.004761905 -0.005295675


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


Re: [R] help: how to change the column name of data.frame

2005-07-19 Thread Prof Brian Ripley
On Tue, 19 Jul 2005, Prof Brian Ripley wrote:

 Just change the names!  E.g.

 names(DF)[c(4,6)] - names(DF)[c(6,4)]

 Strictly a data frame has names, not column names, hence the use the
 names() and names-() functions here.

That answered the subject line.  If you want to exchange the columns (not 
just the data but also the names) you can just use

DF[c(4,6)] - DF[c(6,4)]

Please do try to be more precise as to what you want to do, for example by 
giving an example.

 On Tue, 19 Jul 2005, wu sz wrote:

 I have a data frame with 15 variables, and want to exchange the data
 of 4th column and 6th column. First I append a column in the data
 frame, copy the 4th column data there, then copy the 6th column data
 to 4th column, and copy the appended column data to 6th column, but
 the names of the 4th and 6th column are still unchanged. How can I
 exchange them?

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] Michaelis-menten equation

2005-07-19 Thread Chun-Ying Lee
Dear R users:
   I encountered difficulties in michaelis-menten equation. I found 
that when I use right model definiens, I got wrong Km vlaue, 
and I got right Km value when i use wrong model definiens. 
The value of Vd and Vmax are correct in these two models. 

#-right model definiens
PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
   conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
mm.model - function(time, y, parms) { 
   dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]+y[1]) 
   list(dCpdt)}
Dose-300
modfun - function(time,Vm,Km,Vd) { 
   out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
  rtol=1e-8,atol=1e-8)
  out[,2] } 
objfun - function(par) { 
   out - modfun(PKindex$time,par[1],par[2],par[3]) 
   sum((PKindex$conc-out)^2) } 
fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
print(fit$par)
[1] 10.0390733  0.1341544 34.9891829  #--Km=0.1341544,wrong value--


#-wrong model definiens
#-Km should not divided by Vd--
PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
   conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
mm.model - function(time, y, parms) { 
   dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]/parms[Vd]+y[1]) 
   list(dCpdt)}
Dose-300
modfun - function(time,Vm,Km,Vd) { 
out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
rtol=1e-8,atol=1e-8)
   out[,2] 
} 
objfun - function(par) { 
out - modfun(PKindex$time,par[1],par[2],par[3]) 
sum((PKindex$conc-out)^2)} 
fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
print(fit$par)
[1] 10.038821  4.690267 34.989239  #--Km=4.690267,right value--

What did I do wrong, and how to fix it?
Any suggestions would be greatly appreciated.
Thanks in advance!!

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

Re: [R] Michaelis-menten equation

2005-07-19 Thread Peter Dalgaard
Chun-Ying Lee [EMAIL PROTECTED] writes:

 Dear R users:
I encountered difficulties in michaelis-menten equation. I found 
 that when I use right model definiens, I got wrong Km vlaue, 
 and I got right Km value when i use wrong model definiens. 
 The value of Vd and Vmax are correct in these two models. 

How do you know what the correct value is? Are you sure that the other
values are right?

I'm a bit rusty on MM, but are you sure your right model is right?
Try doing a dimensional analysis on the ODE. I kind of suspect that
Vd is entering in the wrong way. Since you're dealing in
concentrations, should it enter at all (except via the conc. at time
0, of course)?

Not knowing the context, I can't be quite sure, but generally, I'd
expect Vm*Km/(Km+y) to be the reaction rate, so that Vm is the maximum
rate, attained when y is zero and Km is the conc. at half-maximum
rate. This doesn't look quit like what you have. 
 
 #-right model definiens
 PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
 mm.model - function(time, y, parms) { 
dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]+y[1]) 
list(dCpdt)}
 Dose-300
 modfun - function(time,Vm,Km,Vd) { 
out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
   rtol=1e-8,atol=1e-8)
   out[,2] } 
 objfun - function(par) { 
out - modfun(PKindex$time,par[1],par[2],par[3]) 
sum((PKindex$conc-out)^2) } 
 fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
 print(fit$par)
 [1] 10.0390733  0.1341544 34.9891829  #--Km=0.1341544,wrong value--
 
 
 #-wrong model definiens
 #-Km should not divided by Vd--
 PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
 mm.model - function(time, y, parms) { 
dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]/parms[Vd]+y[1]) 
list(dCpdt)}
 Dose-300
 modfun - function(time,Vm,Km,Vd) { 
 out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
 rtol=1e-8,atol=1e-8)
out[,2] 
 } 
 objfun - function(par) { 
 out - modfun(PKindex$time,par[1],par[2],par[3]) 
 sum((PKindex$conc-out)^2)} 
 fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
 print(fit$par)
 [1] 10.038821  4.690267 34.989239  #--Km=4.690267,right value--
 
 What did I do wrong, and how to fix it?
 Any suggestions would be greatly appreciated.
 Thanks in advance!!
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
   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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

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


Re: [R] help: how to change the column name of data.frame

2005-07-19 Thread wu sz
Thanks a lot! But  DF[c(4,6)] - DF[c(6,4)]   seems to just exchange
the data, not the names. If exchanging the columns(both data and
names) needs two steps?

DF[c(4,6)] - DF[c(6,4)]
names(DF)[c(4,6)] - names(DF)[c(6,4)]

Shengzhe

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


Re: [R] help: how to change the column name of data.frame

2005-07-19 Thread Prof Brian Ripley
On Tue, 19 Jul 2005, wu sz wrote:

 Thanks a lot! But  DF[c(4,6)] - DF[c(6,4)]   seems to just exchange
 the data, not the names. If exchanging the columns(both data and
 names) needs two steps?

 DF[c(4,6)] - DF[c(6,4)]
 names(DF)[c(4,6)] - names(DF)[c(6,4)]

Yes, it does.  You did however say in your first message that you wanted 
to exchange the data, despite the subject line.

As I asked before *PLEASE* do try to be precise in what you want to do.

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Michaelis-menten equation

2005-07-19 Thread joerg van den hoff
Chun-Ying Lee wrote:
 Dear R users:
I encountered difficulties in michaelis-menten equation. I found 
 that when I use right model definiens, I got wrong Km vlaue, 
 and I got right Km value when i use wrong model definiens. 
 The value of Vd and Vmax are correct in these two models. 
 
 #-right model definiens
 PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
 mm.model - function(time, y, parms) { 
dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]+y[1]) 
list(dCpdt)}
 Dose-300
 modfun - function(time,Vm,Km,Vd) { 
out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
   rtol=1e-8,atol=1e-8)
   out[,2] } 
 objfun - function(par) { 
out - modfun(PKindex$time,par[1],par[2],par[3]) 
sum((PKindex$conc-out)^2) } 
 fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
 print(fit$par)
 [1] 10.0390733  0.1341544 34.9891829  #--Km=0.1341544,wrong value--
 
 
 #-wrong model definiens
 #-Km should not divided by Vd--
 PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
 mm.model - function(time, y, parms) { 
dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]/parms[Vd]+y[1]) 
list(dCpdt)}
 Dose-300
 modfun - function(time,Vm,Km,Vd) { 
 out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
 rtol=1e-8,atol=1e-8)
out[,2] 
 } 
 objfun - function(par) { 
 out - modfun(PKindex$time,par[1],par[2],par[3]) 
 sum((PKindex$conc-out)^2)} 
 fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
 print(fit$par)
 [1] 10.038821  4.690267 34.989239  #--Km=4.690267,right value--
 
 What did I do wrong, and how to fix it?
 Any suggestions would be greatly appreciated.
 Thanks in advance!!
 
 
 

it is not clear to me what you are trying to do:
you seem to have a time-concentration-curve in PKindex and you seem to
set up a derivative of this time dependency according
to some model in dCpdt. AFAIKS this scenario is  not directly related to
the situation described by the Michaelis-Menten-Equation which relates
some input concentration with some product concentration. If Vm and
Km are meant to be the canonical symbols,
what is Vd, a volume of distribution? it is impossible to see (at least
for me) what exactly you want to achieve.

(and in any case, I would prefer nls for a least squares fit instead
of 'optim').

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

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

Re: [R] colnames

2005-07-19 Thread Adaikalavan Ramasamy
First, your problem could be boiled down to the following example. See
how the colnames of the two outputs vary.

df - cbind.data.frame( 100=1:2, 200=3:4 )
df/df
  X100 X200
111
211

m  - as.matrix( df )   # coerce to matrix class
m/m
  100 200
1   1   1
2   1   1

It appears that whenever R has to create a new dataframe automatically,
it tries to get nice colnames. See help(data.frame). I am not exactly
sure why this behaviour is different when creating a matrix. But I do
not think this is a major problem for most people. If you coerce your
input to matrix, the problem goes away.


Next, note the following points :
 a) mat[ 1:3, 1:ncol(mat) ] is equivalent to simply mat[ 1:3,  ]. 
 b) mat[ 2:nrow(mat), ] is equivalent to simply mat[ -1,  ]
See help(subset) for more information.

Using the points above, we can simplify your function as 

 p.RIs2Returns - function (mat){

   mat - as.matrix(mat)
   x - mat[ -nrow(mat), ]
   y - mat[ -1, ]
  
   return( y/x -1 )
 }

If your data contains only numerical data, it is probably good idea to
work with matrices as matrix operations are faster.


Finally, we can shorten your function. You can use the diff (which works
column-wise if input is a matrix) and apply function if you know that 

y/x  =  exp(log(y/x))  =  exp( log(y) - log(x) )

which could be coded in R as

exp( diff( log(r1) ) )

and then subtract 1 from above to get your returns.

Regards, Adai



On Tue, 2005-07-19 at 09:17 +0100, Gilbert Wu wrote:
 Hi Adai,
 
 Many Thanks for the examples.
 
 I work for a financial institution. We are exploring R as a tool to implement 
 our portfolio optimization strategies. Hence, R is still a new language to us.
 
 The script I wrote tried to make a returns matrix from the daily return 
 indices extracted from a SQL database. Please find below the output that 
 produces the 'X' prefix in the colnames. The reason to preserve the column 
 names is that they are stock identifiers which are to be used by other sub 
 systems rather than R.
 
 I would welcome any suggestion to improve the script.
 
 
 Regards,
 
 Gilbert
 
  p.RIs2Returns -
 + function (RIm)
 + {
 + x-RIm[1:(nrow(RIm)-1), 1:ncol(RIm)]
 + y-RIm[2:nrow(RIm), 1:ncol(RIm)]
 + RReturns - (y/x -1)
 + RReturns
 + }
  
  
  channel-odbcConnect(ourSQLDB)
  result-sqlQuery(channel,paste(select * from equityRIs;))
  odbcClose(channel)
  result
stockidsdate  dbPrice
 1   899188 20050713  7.59500
 2   899188 20050714  7.60500
 3   899188 20050715  7.48000
 4   899188 20050718  7.41500
 5   902232 20050713 10.97000
 6   902232 20050714 10.94000
 7   902232 20050715 10.99000
 8   902232 20050718 11.05000
 9   901714 20050713 17.96999
 10  901714 20050714 18.00999
 11  901714 20050715 17.64999
 12  901714 20050718 17.64000
 13  28176U 20050713  5.19250
 14  28176U 20050714  5.25000
 15  28176U 20050715  5.25000
 16  28176U 20050718  5.22500
 17  15322M 20050713 11.44000
 18  15322M 20050714 11.5
 19  15322M 20050715 11.33000
 20  15322M 20050718 11.27000
  r1-reshape(result, timevar=stockid, idvar=sdate, direction=wide)
  r1
  sdate dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
 dbPrice.15322M
 1 20050713  7.595  10.97   17.96999 5.1925
   11.44
 2 20050714  7.605  10.94   18.00999 5.2500
   11.50
 3 20050715  7.480  10.99   17.64999 5.2500
   11.33
 4 20050718  7.415  11.05   17.64000 5.2250
   11.27
  #Set sdate as the rownames
  rownames(r1) -as.character(r1[1:nrow(r1),1:1])
  #Get rid of the first column
  r1 - r1[1:nrow(r1),2:ncol(r1)]
  r1
  dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
 dbPrice.15322M
 20050713  7.595  10.97   17.96999 5.1925  
 11.44
 20050714  7.605  10.94   18.00999 5.2500  
 11.50
 20050715  7.480  10.99   17.64999 5.2500  
 11.33
 20050718  7.415  11.05   17.64000 5.2250  
 11.27
  colnames(r1) - as.character(sub([[:alnum:]]*\\.,, colnames(r1)))
  r1
  899188 902232   901714 28176U 15322M
 20050713  7.595  10.97 17.96999 5.1925  11.44
 20050714  7.605  10.94 18.00999 5.2500  11.50
 20050715  7.480  10.99 17.64999 5.2500  11.33
 20050718  7.415  11.05 17.64000 5.2250  11.27
  RRs-p.RIs2Returns(r1)
  RRs
   X899188  X902232  X901714  X28176U  X15322M
 20050714  0.001316656 -0.002734731  0.002225933  0.011073664  0.005244755
 20050715 -0.016436555  0.004570384 -0.019988906  0.0 -0.014782609
 20050718 -0.008689840  0.005459509 -0.000566006 -0.004761905 -0.005295675
  


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


[R] Predict

2005-07-19 Thread Matthias Eggenberger
When I callculate a linear model, then I can compute via confint the
confidencial intervals. the interval level can be chosen. as result, I get
the parameter of the model according to the interval level. 

On the other hand, I can compute the prediction-values for my model as well
with predict(object, type=c(response) etc.). Here I have also the
possibility to chose a level for the confidential intervals. the output are
the calculatet values for the fit, the lower and upper level. 

the problem now is, that when I calculate the values through the linear
model function with the parameter values I get from confint() an I compare
them with the values I get from predict() these values differ extremely. Why
is that so? Does the command predict() calculate the values through an other
routine? That means the command predict() doesn't use the same parameters to
calculate the prediction-values than the ones given by confint()?

Greetings Matthias

-- 
GMX DSL = Maximale Leistung zum minimalen Preis!

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


Re: [R] Survival dummy variables and some questions

2005-07-19 Thread Adaikalavan Ramasamy
Stephen,

This has nothing to do with your R but to do with your email settings. 

Are you sure you are sending mails in plain text ? Your email on the
R-help mailing archive (see link below) appears to be unreadable 
https://stat.ethz.ch/pipermail/r-help/2005-July/074210.html

Please try to use plain text. (See http://expita.com/nomime.html).

Regards, Adai



On Tue, 2005-07-19 at 08:59 +0200, Stephen wrote:
 Many thanks I follow you what you say You can request predicted
 values at any sequence of ages - I guess there are plenty of postings on
 how to do that  Regards, Stephen - Original Message - From:
 Frank E Harrell Jr To: Stephen Cc: Prof Brian Ripley ; Sent:
 Monday, July 18, 2005 6:13 PM Subject: Re: [R] Survival dummy variables
 and some questions  Stephen wrote:  Hi 1. Right perhaps this should
 clarify. I would like to extract  coefficeints for different levels of
 the IVs (covariate). So for  instance, age of onset I would want
 Hazards etc for every 5 years and so  on... The approach I took was to
 categorize the variables (e.g., age of  onset) and then turn theu
 resultant categorical variable into a factor as  opposed to a
 variable... that is when the problems began An  alternative
 approach to pulling out different values at different levels  of the
 variable is what I seek. 2. I looked for the link, but can't   Your
 needs don't require categorization. You can request predicted values 
 at any sequence of ages. If you want hazard ratios you can take 
 differences in predicted log hazards and antilog them.   Frank   
 --  Frank E Harrell Jr Professor and Chair School of Medicine 
 Department of Biostatistics Vanderbilt University  
 
  ??  
 http://mail.nana.co.il
 
   [[alternative HTML version deleted]]
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


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


[R] R: expression

2005-07-19 Thread Clark Allan
hi all

i am having a problem with the expression/paste command

say we estimate a variable, named PHI

it contains the value of say 2

and we want to display this value as  hat(phi) = PHI onto a graphic

i.e.  hat(phi)=2 

how does one do this?

i've tried the following:

1.  legend(-5,.3,expression(hat(phi)*=*PHI))

2.  legend(-5,.3,paste(expression(phi),=,PHI))

but they do not work.

any help?

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

Re: [R] R: expression

2005-07-19 Thread Uwe Ligges
Clark Allan wrote:

 hi all
 
 i am having a problem with the expression/paste command
 
 say we estimate a variable, named PHI
 
 it contains the value of say 2
 
 and we want to display this value as  hat(phi) = PHI onto a graphic
 
 i.e.hat(phi)=2 
 
 how does one do this?
 
 i've tried the following:
 
 1.legend(-5,.3,expression(hat(phi)*=*PHI))
 
 2.legend(-5,.3,paste(expression(phi),=,PHI))


See ?plotmath or the Help Desk Article Automation of Mathematical 
Annotation in Plots in R News 2 (3), 32-34.

 legend(-5, .3, substitute(hat(phi) == PHI, list(PHI = PHI)))

Uwe Ligges


 but they do not work.
 
 any help?
 
 /
 allan
 
 
 
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

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


Re: [R] R: expression

2005-07-19 Thread Adaikalavan Ramasamy
Something like this :

 x - 0.5
 plot( 1:10, main=substitute( hat(Phi) ~ = ~ x, list(x=x) ) )

Also see http://tolstoy.newcastle.edu.au/R/help/04/09/3371.html

Regards, Adai



On Tue, 2005-07-19 at 13:35 +0200, Clark Allan wrote:
 hi all
 
 i am having a problem with the expression/paste command
 
 say we estimate a variable, named PHI
 
 it contains the value of say 2
 
 and we want to display this value as  hat(phi) = PHI onto a graphic
 
 i.e.hat(phi)=2 
 
 how does one do this?
 
 i've tried the following:
 
 1.legend(-5,.3,expression(hat(phi)*=*PHI))
 
 2.legend(-5,.3,paste(expression(phi),=,PHI))
 
 but they do not work.
 
 any help?
 
 /
 allan
 __ R-help@stat.math.ethz.ch 
 mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the 
 posting guide! http://www.R-project.org/posting-guide.html

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


Re: [R] dataframes of unequal size

2005-07-19 Thread Liaw, Andy
Seems like no one had responded to this one yet, so I'll take a stab:

## Generate some bogus data:
set.seed(45)
dat - cbind(expand.grid(LETTERS[1:2], 1:3), round(runif(6), 2))
names(dat) - c(state, psu, weight)
dat2 - data.frame(state=sample(c(A, B), 100, replace=TRUE),
   psu=sample(3, 100, replace=TRUE),
   weight=rep(0, 100))

## The actual work:
split(dat2$weight, interaction(dat2$state, dat2$psu)) -
split(dat$weight, interaction(dat$state, dat$psu))

This, I think, will only work correctly if all state/psu combinations in
your C are also present in C1.  If not, you can just augment C1 to
include them.

HTH,
Andy


 From: Renuka Sane
 
 I have two dataframes C and C1. Each has three columns viz. state, psu
 and weight. The dataframes are of unequal size i.e. C1 could be
 2/25/50 rows and C has 42000 rows.  C1 is the master table i.e.
 C1$state, C1$psu and C1$weight are never the same. ThisA. P., Urban, 0
 is not so for C.
 
 For example
 C
 state, psu,weight
 A. P., Urban, 0
 Mah., Rural, 0
 W.B., Rural,0
 Ass., Rural,0
 M. P., Urban,0
 A. P., Urban, 0
 ...
 
 C1
 state, psu, weight
 A. P., Urban, 1.3
 A. P., Rural, 1.2
 M. P., Urban, 0.8
 ..
 
 For every row of C, I want to check if C$state==C1$state and
 C$psu==C1$psu. If it is, I want C$weight - C1$weight, else C$weight
 should be zero.
 
 I am doing the following
 for( i in 1:length(C$weight)) {
  C$w[C$state[i]==C1$state  C$psu[i]==C1$psu] - C1$w[C$state[i] ==
 C1$state  C$psu[i] == C1$psu]
 }
 
 This gives me the correct replacements for the number of rows in C1
 and then just repeats the same weights for the remaning rows in C.
 
 Can someone point out the error in what I am doing or show the correct
 way of doing this?
 
 Thanks,
 Renuka
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! 
 http://www.R-project.org/posting-guide.html
 
 


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


[R] data mining

2005-07-19 Thread secretario academico FACEA

Dear all,
I'm looking for some material on data mining with R. I have something 
from Luis Torgo but I'd like to see something else.

If anybody could help me I'll be thankful
Adrián
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] reading data

2005-07-19 Thread secretario academico FACEA

Dear all,
How can I read data from posgresql?
Thanks
Adrián
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] lme outer

2005-07-19 Thread yann clough
Dear R users

a question about outer explanatory variables in lme:

I have measured size of a population of insects in fields.
These fields were spread out over a large region.
The fields are grouped (spatially) in pairs: one with fertiliser high, 
the other one low.
I want to test effect of mean temperature and fertiliser on popsize.
Meantemp was measured for each field, but measurements are correlated 
within pairs,
and this should be taken into account to avoid pseudoreplication (in 
other words, I d like meantemp to be considered an outer variable).

Do I need to replace the temperature measurements by the means for each 
pair?
Or can I leave in the measurements per field pair?

this is my data and model (with meantemp values for each field):
popsize=c(8,19,13,28,30,29,45,41,21,30,20,32,44,52,65,45)
meantemp=c(10,10.4,11.2,11.4,12,12.25,12.5,12.7,10.1,10.7,11.5,11.3,11.7,12.3,12.9,12.8)
fertiliser=as.factor(rep(c(low,high),each=8))
pair=as.factor(rep(c(1:8),times=2))

model1=lme(popsize~meantemp+fertiliser, random=~1|pair)

I now create a vector with the values of meantemp averaged per pair

meantemp2=tapply(meantemp,pair,mean)
meantemp2=meantemp2[pair]

rerun a model with that explanatory variable:
model2=lme(popsize~meantemp2+fertiliser, random=~1|pair)

summary.lme and anova.lme suggest minute differences in the estimated 
parameters and DF (!) between model1 and model2.
How do I explain these differences, especially in the DF?
Is there a model to prefer?

Sincerely,
Yann
-- 
Yann Clough
Agroecology
Georg-August University
Waldweg 26
D-37073 Goettingen
Tel: 0551/39-2358
email: [EMAIL PROTECTED]
www: http://wwwuser.gwdg.de/~uaoe/mitarbeiter/y_clough_e.htm

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


Re: [R] data mining

2005-07-19 Thread Achim Zeileis
On Tue, 19 Jul 2005 10:12:48 -0300 secretario academico FACEA wrote:

 Dear all,
 I'm looking for some material on data mining with R. I have something 
 from Luis Torgo but I'd like to see something else.

There have been several discussions on the list, so browsing the
archives will probably bring up some helpful pointers. Also have a look
at the MachineLearning view at
  http://CRAN.R-project.org/src/contrib/Views/MachineLearning.html
Z

 If anybody could help me I'll be thankful
 Adrián


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


[R] R: stats

2005-07-19 Thread Clark Allan
for the stats gurus

does anyone know if there exists a general formula relating the median
of a continuous distribution to its moments. the distribution could be
skewed or symmetric and is definitely not normal.

the reason for asking is since the median of the particular distribution
that i am interested in is difficult (probably impossible) to obtain.
the median depends depends on an incomplete gamma distribution. the
moments however can be obtained fairly easily.

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

Re: [R] R: stats

2005-07-19 Thread Duncan Murdoch
On 7/19/2005 9:28 AM, Clark Allan wrote:
 for the stats gurus
 
 does anyone know if there exists a general formula relating the median
 of a continuous distribution to its moments. the distribution could be
 skewed or symmetric and is definitely not normal.

Not in general, and probably not any more practically useful than 
F^(-1)(0.5).

 the reason for asking is since the median of the particular distribution
 that i am interested in is difficult (probably impossible) to obtain.
 the median depends depends on an incomplete gamma distribution. the
 moments however can be obtained fairly easily.

R can calculate the incomplete gamma function (see ?pgamma), so that's 
not necessarily a stumbling block.

Duncan Murdoch

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


Re: [R] R: expression

2005-07-19 Thread Thomas Lumley
On Tue, 19 Jul 2005, Uwe Ligges wrote:

 Clark Allan wrote:

 hi all
 i am having a problem with the expression/paste command
 say we estimate a variable, named PHI
 it contains the value of say 2
 and we want to display this value as  hat(phi) = PHI onto a graphic
 i.e.   hat(phi)=2 
 how does one do this?



 legend(-5, .3, substitute(hat(phi) == PHI, list(PHI = PHI)))


or
  legend(-5, .3, bquote(hat(phi) == .(PHI)))


-thomas

Thomas Lumley   Assoc. Professor, Biostatistics
[EMAIL PROTECTED]   University of Washington, Seattle

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


[R] R Site Search is back up

2005-07-19 Thread Jonathan Baron
The search site at
http://finzi.psych.upenn.edu/
is back up.  This allows you to search mailing list archives,
R functions from most packages, and documents.

You can also use the R function RSiteSearch() in R itself, which
opens the results in your browser.

I'm sorry the site was down for so long.  I plan to work out a
system that will allow faster restoration if this happens again
(mutual complete backups with a colleague).  Backups are useless
if you don't have a machine to put them on, and it took 5 days to
get one, then one day to restore.

Some details that you don't need to read:

The site was reconstructed from scratch.  This means that links
from one message in the archives to another message in the
archives will not work anymore.  They never worked very well,
since the message numbers have been changed a few times over the
years.  (Threads work.  Just links don't work.)

In my rush to get this back up, I left out a few packages that
would not install.  I plan to install at least the help files
from these in due course.  I'm also not sure about Bioconductor
packages that are not already on CRAN.  I left out all of these.
Let me know if you want them.  Same with Jim Lindsey's packages.

I also left out all the mail from r-sig-geo, which I had been
including, although it seems hardly worth the effort since volume
is so low now on that list.  If you want it back, let me know.

What happened is still not clear.  There were several things
wrong with the computer, still not fully diagnosed.  The symptoms 
are consistent with overheating, but the fan was working fine,
and so has the air conditioning in the building (which has been
working very hard, global warming having finally hit
Philadelphia).  Recovery from the disk was possible but very
expensive.

What I learned from this is that backup systems do no good unless 
you have a computer to use.  And also I could have backed up in a 
way that made it easier for myself to reconstruct, but that was
not the major source of the delay.

Jon
-- 
Jonathan Baron, Professor of Psychology, University of Pennsylvania
Home page: http://www.sas.upenn.edu/~baron

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


Re: [R] colnames

2005-07-19 Thread Prof Brian Ripley
On Tue, 19 Jul 2005, Adaikalavan Ramasamy wrote:

 First, your problem could be boiled down to the following example. See
 how the colnames of the two outputs vary.

 df - cbind.data.frame( 100=1:2, 200=3:4 )
 df/df
  X100 X200
 111
 211

That one is probably unintentional.

 m  - as.matrix( df )   # coerce to matrix class
 m/m
  100 200
 1   1   1
 2   1   1

 It appears that whenever R has to create a new dataframe automatically,
 it tries to get nice colnames. See help(data.frame). I am not exactly
 sure why this behaviour is different when creating a matrix. But I do
 not think this is a major problem for most people. If you coerce your
 input to matrix, the problem goes away.

A data frame is column-oriented, and can be used as a source of variables, 
e.g. by attach() and the data= argument of all the model-fitting 
functions.  That is not the purpose of matrices, but is why data.frames 
are made to have syntactic names for columns by default.

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
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@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] colnames

2005-07-19 Thread Gilbert Wu
Hi Adai,

Thank you very much for your suggestions. Your optimized function would come in 
very handy cause I will need to generate a matrix of size around 2250 * 1000.

Regards,

Gilbert


-Original Message-
From: Adaikalavan Ramasamy [mailto:[EMAIL PROTECTED]
Sent: 19 July 2005 12:20
To: Gilbert Wu
Cc: r-help@stat.math.ethz.ch
Subject: RE: [R] colnames


First, your problem could be boiled down to the following example. See
how the colnames of the two outputs vary.

df - cbind.data.frame( 100=1:2, 200=3:4 )
df/df
  X100 X200
111
211

m  - as.matrix( df )   # coerce to matrix class
m/m
  100 200
1   1   1
2   1   1

It appears that whenever R has to create a new dataframe automatically,
it tries to get nice colnames. See help(data.frame). I am not exactly
sure why this behaviour is different when creating a matrix. But I do
not think this is a major problem for most people. If you coerce your
input to matrix, the problem goes away.


Next, note the following points :
 a) mat[ 1:3, 1:ncol(mat) ] is equivalent to simply mat[ 1:3,  ]. 
 b) mat[ 2:nrow(mat), ] is equivalent to simply mat[ -1,  ]
See help(subset) for more information.

Using the points above, we can simplify your function as 

 p.RIs2Returns - function (mat){

   mat - as.matrix(mat)
   x - mat[ -nrow(mat), ]
   y - mat[ -1, ]
  
   return( y/x -1 )
 }

If your data contains only numerical data, it is probably good idea to
work with matrices as matrix operations are faster.


Finally, we can shorten your function. You can use the diff (which works
column-wise if input is a matrix) and apply function if you know that 

y/x  =  exp(log(y/x))  =  exp( log(y) - log(x) )

which could be coded in R as

exp( diff( log(r1) ) )

and then subtract 1 from above to get your returns.

Regards, Adai



On Tue, 2005-07-19 at 09:17 +0100, Gilbert Wu wrote:
 Hi Adai,
 
 Many Thanks for the examples.
 
 I work for a financial institution. We are exploring R as a tool to implement 
 our portfolio optimization strategies. Hence, R is still a new language to us.
 
 The script I wrote tried to make a returns matrix from the daily return 
 indices extracted from a SQL database. Please find below the output that 
 produces the 'X' prefix in the colnames. The reason to preserve the column 
 names is that they are stock identifiers which are to be used by other sub 
 systems rather than R.
 
 I would welcome any suggestion to improve the script.
 
 
 Regards,
 
 Gilbert
 
  p.RIs2Returns -
 + function (RIm)
 + {
 + x-RIm[1:(nrow(RIm)-1), 1:ncol(RIm)]
 + y-RIm[2:nrow(RIm), 1:ncol(RIm)]
 + RReturns - (y/x -1)
 + RReturns
 + }
  
  
  channel-odbcConnect(ourSQLDB)
  result-sqlQuery(channel,paste(select * from equityRIs;))
  odbcClose(channel)
  result
stockidsdate  dbPrice
 1   899188 20050713  7.59500
 2   899188 20050714  7.60500
 3   899188 20050715  7.48000
 4   899188 20050718  7.41500
 5   902232 20050713 10.97000
 6   902232 20050714 10.94000
 7   902232 20050715 10.99000
 8   902232 20050718 11.05000
 9   901714 20050713 17.96999
 10  901714 20050714 18.00999
 11  901714 20050715 17.64999
 12  901714 20050718 17.64000
 13  28176U 20050713  5.19250
 14  28176U 20050714  5.25000
 15  28176U 20050715  5.25000
 16  28176U 20050718  5.22500
 17  15322M 20050713 11.44000
 18  15322M 20050714 11.5
 19  15322M 20050715 11.33000
 20  15322M 20050718 11.27000
  r1-reshape(result, timevar=stockid, idvar=sdate, direction=wide)
  r1
  sdate dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
 dbPrice.15322M
 1 20050713  7.595  10.97   17.96999 5.1925
   11.44
 2 20050714  7.605  10.94   18.00999 5.2500
   11.50
 3 20050715  7.480  10.99   17.64999 5.2500
   11.33
 4 20050718  7.415  11.05   17.64000 5.2250
   11.27
  #Set sdate as the rownames
  rownames(r1) -as.character(r1[1:nrow(r1),1:1])
  #Get rid of the first column
  r1 - r1[1:nrow(r1),2:ncol(r1)]
  r1
  dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
 dbPrice.15322M
 20050713  7.595  10.97   17.96999 5.1925  
 11.44
 20050714  7.605  10.94   18.00999 5.2500  
 11.50
 20050715  7.480  10.99   17.64999 5.2500  
 11.33
 20050718  7.415  11.05   17.64000 5.2250  
 11.27
  colnames(r1) - as.character(sub([[:alnum:]]*\\.,, colnames(r1)))
  r1
  899188 902232   901714 28176U 15322M
 20050713  7.595  10.97 17.96999 5.1925  11.44
 20050714  7.605  10.94 18.00999 5.2500  11.50
 20050715  7.480  10.99 17.64999 5.2500  11.33
 20050718  7.415  11.05 17.64000 5.2250  11.27
  RRs-p.RIs2Returns(r1)
  RRs
   X899188  X902232  X901714  X28176U  X15322M
 20050714  0.001316656 -0.002734731  0.002225933  0.011073664  0.005244755
 

Re: [R] Survival dummy variables and some questions

2005-07-19 Thread Stephen
Thanks For pointing that out. S - Original Message - From:
Adaikalavan Ramasamy To: Stephen Cc: Frank E Harrell Jr ; Prof
Brian Ripley ; Sent: Tuesday, July 19, 2005 1:30 PM Subject: Re: [R]
Survival dummy variables and some questions  Stephen,   This has
nothing to do with your R but to do with your email settings.   Are
you sure you are sending mails in plain text ? Your email on the 
R-help mailing archive (see link below) appears to be unreadable 
https://stat.ethz.ch/pipermail/r-help/2005-July/074210.html   Please
try to use plain text. (See http://expita.com/nomime.html).   Regards,
Adai On Tue, 2005-07-19 at 08:59 +0200, Stephen wrote:  Many
thanks I follow you what you say You can request predicted  values
at any sequence of ages - I guess there are plenty of postings on  how
to do that  Regards, Stephen - Original Message - From: 
Frank E Harrell Jr To: Stephen Cc: Prof Brian Ripley ; Sent: 
Monday, July 18, 2005 6:13 PM Subject: Re: [R] Survival dummy variables
 and some questions  Stephen wrote:  Hi 1. Right perhaps this
should  clarify. I would like to extract  coefficeints for different
levels of  the IVs (covariate). So for  instance, age of onset I
would want  Hazards etc for every 5 years and so  on... The approach
I took was to  categorize the variables (e.g., age of  onset) and
then turn theu  resultant categorical variable into a factor as 
opposed to a  variable... that is when the problems began An 
alternative  approach to pulling out different values at different
levels  of the  variable is what I seek. 2. I looked for the link,
but can't   Your  needs don't require categorization. You can
request predicted values   at any sequence of ages. If you want
hazard ratios you can take   differences in predicted log hazards and
antilog them.   Frank --  Frank E Harrell Jr Professor and
Chair School of Medicine   Department of Biostatistics Vanderbilt
University ??    http://mail.nana.co.il  
[[alternative HTML version deleted]]  
__ 
R-help@stat.math.ethz.ch mailing list 
https://stat.ethz.ch/mailman/listinfo/r-help  PLEASE do read the
posting guide!  http://www.R-project.org/posting-guide.html

 ??  
http://mail.nana.co.il

[[alternative HTML version deleted]]

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


[R] extracting row means from a list

2005-07-19 Thread Andy Bunn
Hello: I'm reading in a series of text files (100 files that are each 2000
rows by 6 columns). I wish to combine the columns (6) of each file (100) and
get the row mean. I'd like to end up with a data.frame of 2000 rows by 6
columns.

foo - list()
for(i in 1:10){
 # The real data are read in from a series of numbered text files
 foo[[i]] - data.frame(x1 = rnorm(100), x2 = rnorm(100), x3 =
rnorm(100),
x4 = rnorm(100), x5 = rnorm(100), x6 =
rnorm(100))

}

str(foo)
# by hand
mean.x1 -
rowMeans(cbind(foo[[1]][,1],foo[[2]][,1],foo[[3]][,1],foo[[4]][,1],foo[[5]][
,1]),
  foo[[6]][,1],foo[[7]][,1],foo[[8]][,1],foo[[9]][,1
],foo[[10]][,1]))
mean.x2 -
rowMeans(cbind(foo[[1]][,2],foo[[2]][,2],foo[[3]][,2],foo[[4]][,2],foo[[5]][
,2]),
  foo[[6]][,2],foo[[7]][,2],foo[[8]][,2],foo[[9]][,2
],foo[[10]][,2]))
# and so on to column 6
mean.x6 -
rowMeans(cbind(foo[[1]][,6],foo[[2]][,6],foo[[3]][,6],foo[[4]][,6],foo[[5]][
,6]),
  foo[[6]][,6],foo[[7]][,6],foo[[8]][,6],foo[[9]][,6
],foo[[10]][,6]))


I've implemented this with nested loops that create temporary variables and
calc the mean, but the approach is clunky. E.g.,

# nested loops
for(i in 1:ncol(foo[[1]])){
  for(j in 1:length(foo)){
# etc ...
  }
}

Is there a way to build a better mouse trap?

TIA, Andy

Thanks, Andy

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


[R] using argument names (of indeterminate number) within a function

2005-07-19 Thread Dirk Enzmann
Although I tried to find an answer in the manuals and archives, I cannot 
solve this (please excuse that my English and/or R programming skills 
are not good enough to state my problem more clearly):

I want to write a function with an indeterminate (not pre-defined) 
number of arguments and think that I should use the ... construct and 
the match.call() function. The goal is to write a function that (among 
other things) uses cbind() to combine a not pre-defined number of 
vectors specified in the function call. For example, if my vectors are 
x1, x2, x3, ... xn, within the function I want to use cbind(x1, x2) or 
cbind(x1, x3, x5) or ... depending on the vector names I use in the 
funcion call. Additionally, the function has other arguments.

In the archives I found the following thread (followed by Marc Schwartz)

http://finzi.psych.upenn.edu/R/Rhelp02a/archive/15186.html
[R] returning argument names from Peter Dalgaard BSA on 2003-04-10 (stdin)

that seems to contain the solution to my problem, but I am stuck because 
  sapply(match.call()[-1], deparse) gives me a vector of strings and I 
don't know how to use the names in this vector in the cbind() function.

Up to now my (clearly deficit) function looks like:

test - function(..., mvalid=1)
{
   args = sapply(match.call()[-1], deparse)
# and here, I don't know how the vector names in args
# can be used in the cbind() function to follow:
#
# temp - cbind( ???
   if (mvalid  1)
   {
#  here it goes on
   }
}

Ultimately, I want that the function can be called like
test(x1,x2,mvalid=1)
or
test(x1,x3,x5,mavlid=2)
and that within the function
cbind(x1,x2)
or cbind(x1,x3,x5)
will be used.

Can someone give and explain an example / a solution on how to proceed?

*
Dr. Dirk Enzmann
Institute of Criminal Sciences
Dept. of Criminology
Edmund-Siemers-Allee 1
D-20146 Hamburg
Germany

phone: +49-040-42838.7498 (office)
+49-040-42838.4591 (Billon)
fax:   +49-040-42838.2344
email: [EMAIL PROTECTED]
www: 
http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Enzmann.html

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


Re: [R] colnames

2005-07-19 Thread Gilbert Wu
Hi Adai,

When I tried the optimized routine, I got the following error message:

r1
 899188 902232   901714 28176U 15322M
20050713  7.595  10.97 17.96999 5.1925  11.44
20050714  7.605  10.94 18.00999 5.2500  11.50
20050715  7.480  10.99 17.64999 5.2500  11.33
20050718  7.415  11.05 17.64000 5.2250  11.27
 exp(diff(log(r1))) -1
Error in r[i1] - r[-length(r):-(length(r) - lag + 1)] : 
non-numeric argument to binary operator


Any idea?

Many Thanks.

Gilbert
-Original Message-
From: Adaikalavan Ramasamy [mailto:[EMAIL PROTECTED]
Sent: 19 July 2005 12:20
To: Gilbert Wu
Cc: r-help@stat.math.ethz.ch
Subject: RE: [R] colnames


First, your problem could be boiled down to the following example. See
how the colnames of the two outputs vary.

df - cbind.data.frame( 100=1:2, 200=3:4 )
df/df
  X100 X200
111
211

m  - as.matrix( df )   # coerce to matrix class
m/m
  100 200
1   1   1
2   1   1

It appears that whenever R has to create a new dataframe automatically,
it tries to get nice colnames. See help(data.frame). I am not exactly
sure why this behaviour is different when creating a matrix. But I do
not think this is a major problem for most people. If you coerce your
input to matrix, the problem goes away.


Next, note the following points :
 a) mat[ 1:3, 1:ncol(mat) ] is equivalent to simply mat[ 1:3,  ]. 
 b) mat[ 2:nrow(mat), ] is equivalent to simply mat[ -1,  ]
See help(subset) for more information.

Using the points above, we can simplify your function as 

 p.RIs2Returns - function (mat){

   mat - as.matrix(mat)
   x - mat[ -nrow(mat), ]
   y - mat[ -1, ]
  
   return( y/x -1 )
 }

If your data contains only numerical data, it is probably good idea to
work with matrices as matrix operations are faster.


Finally, we can shorten your function. You can use the diff (which works
column-wise if input is a matrix) and apply function if you know that 

y/x  =  exp(log(y/x))  =  exp( log(y) - log(x) )

which could be coded in R as

exp( diff( log(r1) ) )

and then subtract 1 from above to get your returns.

Regards, Adai



On Tue, 2005-07-19 at 09:17 +0100, Gilbert Wu wrote:
 Hi Adai,
 
 Many Thanks for the examples.
 
 I work for a financial institution. We are exploring R as a tool to implement 
 our portfolio optimization strategies. Hence, R is still a new language to us.
 
 The script I wrote tried to make a returns matrix from the daily return 
 indices extracted from a SQL database. Please find below the output that 
 produces the 'X' prefix in the colnames. The reason to preserve the column 
 names is that they are stock identifiers which are to be used by other sub 
 systems rather than R.
 
 I would welcome any suggestion to improve the script.
 
 
 Regards,
 
 Gilbert
 
  p.RIs2Returns -
 + function (RIm)
 + {
 + x-RIm[1:(nrow(RIm)-1), 1:ncol(RIm)]
 + y-RIm[2:nrow(RIm), 1:ncol(RIm)]
 + RReturns - (y/x -1)
 + RReturns
 + }
  
  
  channel-odbcConnect(ourSQLDB)
  result-sqlQuery(channel,paste(select * from equityRIs;))
  odbcClose(channel)
  result
stockidsdate  dbPrice
 1   899188 20050713  7.59500
 2   899188 20050714  7.60500
 3   899188 20050715  7.48000
 4   899188 20050718  7.41500
 5   902232 20050713 10.97000
 6   902232 20050714 10.94000
 7   902232 20050715 10.99000
 8   902232 20050718 11.05000
 9   901714 20050713 17.96999
 10  901714 20050714 18.00999
 11  901714 20050715 17.64999
 12  901714 20050718 17.64000
 13  28176U 20050713  5.19250
 14  28176U 20050714  5.25000
 15  28176U 20050715  5.25000
 16  28176U 20050718  5.22500
 17  15322M 20050713 11.44000
 18  15322M 20050714 11.5
 19  15322M 20050715 11.33000
 20  15322M 20050718 11.27000
  r1-reshape(result, timevar=stockid, idvar=sdate, direction=wide)
  r1
  sdate dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
 dbPrice.15322M
 1 20050713  7.595  10.97   17.96999 5.1925
   11.44
 2 20050714  7.605  10.94   18.00999 5.2500
   11.50
 3 20050715  7.480  10.99   17.64999 5.2500
   11.33
 4 20050718  7.415  11.05   17.64000 5.2250
   11.27
  #Set sdate as the rownames
  rownames(r1) -as.character(r1[1:nrow(r1),1:1])
  #Get rid of the first column
  r1 - r1[1:nrow(r1),2:ncol(r1)]
  r1
  dbPrice.899188 dbPrice.902232 dbPrice.901714 dbPrice.28176U 
 dbPrice.15322M
 20050713  7.595  10.97   17.96999 5.1925  
 11.44
 20050714  7.605  10.94   18.00999 5.2500  
 11.50
 20050715  7.480  10.99   17.64999 5.2500  
 11.33
 20050718  7.415  11.05   17.64000 5.2250  
 11.27
  colnames(r1) - as.character(sub([[:alnum:]]*\\.,, colnames(r1)))
  r1
  899188 902232   901714 28176U 15322M
 20050713  7.595  10.97 17.96999 5.1925  11.44
 20050714  7.605  10.94 18.00999 5.2500  

Re: [R] extracting row means from a list

2005-07-19 Thread Sundar Dorai-Raj


Andy Bunn wrote:
 Hello: I'm reading in a series of text files (100 files that are each 2000
 rows by 6 columns). I wish to combine the columns (6) of each file (100) and
 get the row mean. I'd like to end up with a data.frame of 2000 rows by 6
 columns.
 
 foo - list()
 for(i in 1:10){
  # The real data are read in from a series of numbered text files
  foo[[i]] - data.frame(x1 = rnorm(100), x2 = rnorm(100), x3 =
 rnorm(100),
 x4 = rnorm(100), x5 = rnorm(100), x6 =
 rnorm(100))
 
 }
 
 str(foo)
 # by hand
 mean.x1 -
 rowMeans(cbind(foo[[1]][,1],foo[[2]][,1],foo[[3]][,1],foo[[4]][,1],foo[[5]][
 ,1]),
   foo[[6]][,1],foo[[7]][,1],foo[[8]][,1],foo[[9]][,1
 ],foo[[10]][,1]))
 mean.x2 -
 rowMeans(cbind(foo[[1]][,2],foo[[2]][,2],foo[[3]][,2],foo[[4]][,2],foo[[5]][
 ,2]),
   foo[[6]][,2],foo[[7]][,2],foo[[8]][,2],foo[[9]][,2
 ],foo[[10]][,2]))
 # and so on to column 6
 mean.x6 -
 rowMeans(cbind(foo[[1]][,6],foo[[2]][,6],foo[[3]][,6],foo[[4]][,6],foo[[5]][
 ,6]),
   foo[[6]][,6],foo[[7]][,6],foo[[8]][,6],foo[[9]][,6
 ],foo[[10]][,6]))
 
 
 I've implemented this with nested loops that create temporary variables and
 calc the mean, but the approach is clunky. E.g.,
 
 # nested loops
 for(i in 1:ncol(foo[[1]])){
   for(j in 1:length(foo)){
 # etc ...
   }
 }
 
 Is there a way to build a better mouse trap?
 
 TIA, Andy


I don't know of a way of getting around at least one for loop, but the 
following might be want you need:

r - matrix(, NROW(foo[[1]]), length(foo))
for(i in 1:NCOL(foo[[1]]))
   r[, i] - rowMeans(do.call(cbind, lapply(foo, [, i)))
dim(r)

This assumes each element of foo has identical dimensions. Otherwise 
you'll get an error.

HTH,

--sundar

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


Re: [R] using argument names (of indeterminate number) within a function

2005-07-19 Thread Marc Schwartz (via MN)
On Tue, 2005-07-19 at 17:29 +0200, Dirk Enzmann wrote:
 Although I tried to find an answer in the manuals and archives, I cannot 
 solve this (please excuse that my English and/or R programming skills 
 are not good enough to state my problem more clearly):
 
 I want to write a function with an indeterminate (not pre-defined) 
 number of arguments and think that I should use the ... construct and 
 the match.call() function. The goal is to write a function that (among 
 other things) uses cbind() to combine a not pre-defined number of 
 vectors specified in the function call. For example, if my vectors are 
 x1, x2, x3, ... xn, within the function I want to use cbind(x1, x2) or 
 cbind(x1, x3, x5) or ... depending on the vector names I use in the 
 funcion call. Additionally, the function has other arguments.
 
 In the archives I found the following thread (followed by Marc Schwartz)
 
 http://finzi.psych.upenn.edu/R/Rhelp02a/archive/15186.html
 [R] returning argument names from Peter Dalgaard BSA on 2003-04-10 (stdin)
 
 that seems to contain the solution to my problem, but I am stuck because 
   sapply(match.call()[-1], deparse) gives me a vector of strings and I 
 don't know how to use the names in this vector in the cbind() function.
 
 Up to now my (clearly deficit) function looks like:
 
 test - function(..., mvalid=1)
 {
args = sapply(match.call()[-1], deparse)
 # and here, I don't know how the vector names in args
 # can be used in the cbind() function to follow:
 #
 # temp - cbind( ???
if (mvalid  1)
{
 #  here it goes on
}
 }
 
 Ultimately, I want that the function can be called like
 test(x1,x2,mvalid=1)
 or
 test(x1,x3,x5,mavlid=2)
 and that within the function
 cbind(x1,x2)
 or cbind(x1,x3,x5)
 will be used.
 
 Can someone give and explain an example / a solution on how to proceed?

Hi Dirk,

How about this:

my.cbind - function(...)
{
  do.call(cbind, list(...))
}

 a - 1:10
 b - 11:20
 c - 21:30
 d - 31:40

 my.cbind(a, b)
  [,1] [,2]
 [1,]1   11
 [2,]2   12
 [3,]3   13
 [4,]4   14
 [5,]5   15
 [6,]6   16
 [7,]7   17
 [8,]8   18
 [9,]9   19
[10,]   10   20

 my.cbind(b, c, d)
  [,1] [,2] [,3]
 [1,]   11   21   31
 [2,]   12   22   32
 [3,]   13   23   33
 [4,]   14   24   34
 [5,]   15   25   35
 [6,]   16   26   36
 [7,]   17   27   37
 [8,]   18   28   38
 [9,]   19   29   39
[10,]   20   30   40

 my.cbind(a, b, c, d)
  [,1] [,2] [,3] [,4]
 [1,]1   11   21   31
 [2,]2   12   22   32
 [3,]3   13   23   33
 [4,]4   14   24   34
 [5,]5   15   25   35
 [6,]6   16   26   36
 [7,]7   17   27   37
 [8,]8   18   28   38
 [9,]9   19   29   39
[10,]   10   20   30   40


The use of list(...) in the function allows you to use list based
functions such as do.call() or lapply() against the argument objects
directly without having to deparse and re-parse the character names of
the arguments, which is the approach that Peter used in his response in
the thread you referenced. 

The OP in that thread wanted the argument names as character vectors, as
opposed to the argument objects themselves, which is what you need here.

HTH,

Marc Schwartz

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


[R] ROC curve with survival data

2005-07-19 Thread SUBIRANA CACHINERO, ISAAC
Hi everyone, 

I am doing 5 years mortality predictive index score with survival analysis 
using a Cox proportional hazard model where I have a continous predictive 
variable and a right censored response which is the mortality, and the 
individuals were followed a maximum of 7 years.

I'd like to asses the discrimination ability of survival analysis Cox model by 
computing a ROC curve and area under the curve for a fixed time (5 years), 
taking into account that the response is not binary but right censored.

So, is there a function that computes a ROC curve under right censored survival 
data?

Thank you in advance.

 

 

Isaac Subirana: [EMAIL PROTECTED]

Institut Municipal d'Investigació Mèdica (IMIM),

Barcelona (Spain)

 


[[alternative HTML version deleted]]

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

[R] .gct file

2005-07-19 Thread mark salsburg
I have two files to compare, one is a regular txt file that I can read
in no prob.

The other is a .gct file (How do I read in this one?)

I tried a simple

read.table(data.gct, header = T)

How do you suggest reading in this file??

thank you.

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


Re: [R] .gct file

2005-07-19 Thread Duncan Murdoch
On 7/19/2005 12:10 PM, mark salsburg wrote:
 I have two files to compare, one is a regular txt file that I can read
 in no prob.
 
 The other is a .gct file (How do I read in this one?)
 
 I tried a simple
 
 read.table(data.gct, header = T)
 
 How do you suggest reading in this file??
 

.gct is not a standard filename extension.  You need to know what is in 
that file.  Where did you get it?  What program created it?

Chances are the easiest thing to do is to get the program that created 
it to export in a well known format, e.g. .csv.

Duncan Murdoch

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


Re: [R] Predict

2005-07-19 Thread Christoph Buser
Dear Matthias

Can you provide an example to demonstrate what you did? Two
remarks to your email. Maybe that answers already your question.

1) Using predict() you will get the estimated value for each
   observation or for new data. You can reproduce this value by
   using the coefficients from your estimated model (see the
   example below). 
   For the interval you can get a confidence interval for the
   expected value under fixed conditions of the explanatory
   variables or you can obtain a prediction interval for a
   single new observation. The latter is of course wider, since
   you try to catch a single observation and not the expected
   value. 

2) Using confint() you will get the estimated parameters (which
   are random variables, too) and their confidence interval. 
   You can use the estimated values to calculate the predicted
   values.

But you can NOT use the upper values from confint to
estimate the upper values from predict by just putting them into
your regression model. Thats not the way how confidence
intervals are constructed.
(I am not sure if this was your intention. Maybe if you show a
reproducible example you can correct me if you meant something
different) 

## R Code
## Creation of a dataframe 
set.seed(1)
x1 - runif(40)
f1 - rep(c(a, b, c,d), each = 10)
y - 2*x1 + rep(c(0.5, 0.1, -0.6, 1.5), each = 10) + rnorm(40, 0, 2)
dat - data.frame(y = y, f1 = f1, x1 = x1)
## regression model
reg - lm(y~ x1 + f1, data = dat)
summary(reg)

confint(reg)
predict(reg, type=c(response), interval = confidence)

## caluclation of predicted values using the estimated
## coefficients 

## estimated coefficients
co - summary(reg)$coefficients[,Estimate]
## Using the regression model with that coefficients
## for observation 11 
co[(Intercept)] + dat[11,x1]*co[x1] + co[f1b]
## prediction of observation 11 
predict(reg, type=c(response))[11]


Regards,

Christoph 

--
Christoph Buser [EMAIL PROTECTED]
Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology)  8092 Zurich  SWITZERLAND
phone: x-41-44-632-4673 fax: 632-1228
http://stat.ethz.ch/~buser/
--



Matthias Eggenberger writes:
  When I callculate a linear model, then I can compute via confint the
  confidencial intervals. the interval level can be chosen. as result, I get
  the parameter of the model according to the interval level. 
  
  On the other hand, I can compute the prediction-values for my model as well
  with predict(object, type=c(response) etc.). Here I have also the
  possibility to chose a level for the confidential intervals. the output are
  the calculatet values for the fit, the lower and upper level. 
  
  the problem now is, that when I calculate the values through the linear
  model function with the parameter values I get from confint() an I compare
  them with the values I get from predict() these values differ extremely. Why
  is that so? Does the command predict() calculate the values through an other
  routine? That means the command predict() doesn't use the same parameters to
  calculate the prediction-values than the ones given by confint()?
  
  Greetings Matthias
  
  -- 
  GMX DSL = Maximale Leistung zum minimalen Preis!
  
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

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


Re: [R] Padding in lattice plots

2005-07-19 Thread Deepayan Sarkar
On 7/16/05, Federico Gherardini [EMAIL PROTECTED] wrote:
 On Friday 15 July 2005 17:00, Deepayan Sarkar wrote:
  On 7/15/05, Federico Gherardini [EMAIL PROTECTED] wrote:
   On Friday 15 July 2005 14:42, you wrote:
Hi all,
I've used the split argument to print four lattice plots on a single
page. The problem now is that I need to reduce the amount of white
space between the plots. I've read other mails in this list about the
new trellis parameters layout.heights and layout.widhts but I haven't
been able to use them properly. I've tried to input values between 0
and 1 as the padding value (both left and right and top and bottom) but
nothing changed. It seems I can only increase the padding by using
values  1. Any ideas?
   
Thanks in advance for your help
Federico Gherardini
  
   It seems like I've found an answer myself you have to use negative
   values to decrease the padding. I thought it was something like the cex
   parameter which acts like a multiplier
 
  I thought so too.
 
   but this is not the case.
 
  Could you post what you used? There are several different padding
  parameters you need to set to 0, did you change them all?
 
  Deepayan
 
 Hi Deepayan
 This is what I used I don't know if I did everything the proper way but
 at least I got the result I was seeking! :)
 
 trellis.par.set(list(layout.heights = list(top.padding = -1)))
 
 trellis.par.set(list(layout.heights = list(bottom.padding = -1,
 axis.xlab.padding = 1, xlab = -1.2)))
 
 trellis.par.set(list(layout.widths = list(left.padding = -1)))
 
 trellis.par.set(list(layout.widths = list(right.padding = -1,
 ylab.axis.padding = -0.5)))
 
 Do these settings make any sense?

Yes, but there are other padding parameters. For example, at the top, there's 

 $ top.padding  : num 1
 $ main.key.padding : num 1
 $ key.axis.padding : num 1

The total amount of space left at the top (in the default case, with
no main label and no key) is the sum of all three, so just setting one
to 0 wouldn't be enough for what you want. (This is probably not very
simple, but I couldn't think of anything that's simpler yet as
flexible.)

Deepayan

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


[R] initial points for arms in package HI

2005-07-19 Thread Christoph Buser
Dear R-users

I have a problem choosing initial points for the function arms()
in the package HI
I intend to implement a Gibbs sampler and one of my conditional
distributions is nonstandard and not logconcave.
Therefore I'd like to use arms.

But there seem to be a strong influence of the initial point
y.start. To show the effect I constructed a demonstration
example. It is reproducible without further information.  

Please note that my target density is not logconcave.

Thanks for all comments or ideas.

Christoph Buser

## R Code:

library(HI)
## parameter for the distribution
para - 0.1

## logdensity
logDichteGam - function(x, u = para, v = para) {
  -(u*x + v*1/x) - log(x)
}
## density except for the constant
propDichteGam - function(x, u = para, v = para) {
  exp(-(u*x + v*1/x) - log(x))
}
## calculating the constant 
(c - integrate(propDichteGam, 0, 1000, rel.tol = 10^(-12))$value)
## density
DichteGam - function(x, u = para, v = para) {
  exp(-(u*x + v*1/x) - log(x))/c
}

## calculating 1000 values by repeating a call of arms (this would
## be the situation in an Gibbs Sample. Of course in a Gibbs sampler
## the distribution would change. This is only for demonstration
res1 - NULL
for(i in 1:1000)
  res1[i] - arms(runif(1,0,100), logDichteGam, function(x) (x0)(x100), 1)

## Generating a sample of thousand observations with 1 call of arms
res2 - arms(runif(1,0,100), logDichteGam, function(x) (x0)(x100), 1000)

## Plot of the samples 
mult.fig(4)
plot(res1, log = y)
plot(res2, log = y)
hist(res1, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
 ylim = c(0,1))
curve(DichteGam, 0,4, add = TRUE, col = 2)
hist(res2, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
 ylim = c(0,1))
curve(DichteGam, 0,4, add = TRUE, col = 2)


## If we repeat the procedure, using the fix intial value 1,
## the situation is even worse
res3 - NULL
for(i in 1:1000)
  res3[i] - arms(1, logDichteGam, function(x) (x0)(x100), 1)

## Generating a sample of thousand observations with 1 call of arms
res4 - arms(1, logDichteGam, function(x) (x0)(x100), 1000)

## Plot of the samples 
par(mfrow = c(2,2))
plot(res3, log = y)
plot(res4, log = y)
hist(res3, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
 ylim = c(0,1))
curve(DichteGam, 0,4, add = TRUE, col = 2)
hist(res4, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
 ylim = c(0,1))
curve(DichteGam, 0,4, add = TRUE, col = 2)


## If I generate the sample in a for-loop (one by one) I do not
## get the correct density. But this is exactly the situation in 
## my Gibbs Sampler. Therfore I am concerned about the correct 
## application of arms

--
Christoph Buser [EMAIL PROTECTED]
Seminar fuer Statistik, LEO C13
ETH (Federal Inst. Technology)  8092 Zurich  SWITZERLAND
phone: x-41-44-632-4673 fax: 632-1228
http://stat.ethz.ch/~buser/

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


Re: [R] .gct file

2005-07-19 Thread Marc Schwartz (via MN)
On Tue, 2005-07-19 at 12:28 -0400, Duncan Murdoch wrote:
 On 7/19/2005 12:10 PM, mark salsburg wrote:
  I have two files to compare, one is a regular txt file that I can read
  in no prob.
  
  The other is a .gct file (How do I read in this one?)
  
  I tried a simple
  
  read.table(data.gct, header = T)
  
  How do you suggest reading in this file??
  
 
 .gct is not a standard filename extension.  You need to know what is in 
 that file.  Where did you get it?  What program created it?
 
 Chances are the easiest thing to do is to get the program that created 
 it to export in a well known format, e.g. .csv.
 
 Duncan Murdoch

A quick Google search would suggest Gene Cluster Text file:

http://www.broad.mit.edu/cancer/software/genepattern/tutorial/gp_tutorial_fileformats.html#gct

produced by Gene Pattern:

http://www.broad.mit.edu/cancer/software/genepattern/

If correct, I would point Mark to the Bioconductor folks for more
information and assistance:

http://www.bioconductor.org/

HTH,

Marc Schwartz

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


Re: [R] .gct file

2005-07-19 Thread mark salsburg
ok so the gct file looks like this:

#1.2  (version number)
7283 19   (matrix size)
Name Description Values
  ...  ..

How can I tell R to disregard the first two lines and start reading
the 3rd line in this gct file. I would just delete them, but I do not
know how to open a gct. file

thank you

On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
 On 7/19/2005 12:10 PM, mark salsburg wrote:
  I have two files to compare, one is a regular txt file that I can read
  in no prob.
 
  The other is a .gct file (How do I read in this one?)
 
  I tried a simple
 
  read.table(data.gct, header = T)
 
  How do you suggest reading in this file??
 
 
 .gct is not a standard filename extension.  You need to know what is in
 that file.  Where did you get it?  What program created it?
 
 Chances are the easiest thing to do is to get the program that created
 it to export in a well known format, e.g. .csv.
 
 Duncan Murdoch


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


Re: [R] .gct file

2005-07-19 Thread Randy Johnson
If it is a text file ?read.table should provide enough details to read the
file into R. Based on the file format referenced below it shouldn't be too
hard to get at the parts you want.

Randy


On 7/19/05 1:06 PM, Marc Schwartz (via MN) [EMAIL PROTECTED] wrote:

 On Tue, 2005-07-19 at 12:28 -0400, Duncan Murdoch wrote:
 On 7/19/2005 12:10 PM, mark salsburg wrote:
 I have two files to compare, one is a regular txt file that I can read
 in no prob.
 
 The other is a .gct file (How do I read in this one?)
 
 I tried a simple
 
 read.table(data.gct, header = T)
 
 How do you suggest reading in this file??
 
 
 .gct is not a standard filename extension.  You need to know what is in
 that file.  Where did you get it?  What program created it?
 
 Chances are the easiest thing to do is to get the program that created
 it to export in a well known format, e.g. .csv.
 
 Duncan Murdoch
 
 A quick Google search would suggest Gene Cluster Text file:
 
 http://www.broad.mit.edu/cancer/software/genepattern/tutorial/gp_tutorial_file
 formats.html#gct
 
 produced by Gene Pattern:
 
 http://www.broad.mit.edu/cancer/software/genepattern/
 
 If correct, I would point Mark to the Bioconductor folks for more
 information and assistance:
 
 http://www.bioconductor.org/
 
 HTH,
 
 Marc Schwartz
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 

~~
Randy Johnson
Laboratory of Genomic Diversity
NCI-Frederick
Bldg 560, Rm 11-85
Frederick, MD 21702
(301)846-1304
~~

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


Re: [R] .gct file

2005-07-19 Thread Duncan Murdoch
On 7/19/2005 1:16 PM, mark salsburg wrote:
 ok so the gct file looks like this:
 
 #1.2  (version number)
 7283 19   (matrix size)
 Name Description Values
   ...  ..
 
 How can I tell R to disregard the first two lines and start reading
 the 3rd line in this gct file. I would just delete them, but I do not
 know how to open a gct. file

Use skip=2.  See ?read.table.

Duncan Murdoch

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


[R] Library mclust in 64bit compiled R

2005-07-19 Thread Tae-Hoon Chung
Hi, All;

I tried to use library mclust in 64-bit compiled R 2.0.1 but failed.
Installation went smoothly without any warning or error. However, when I
tried to use them with the following simple code, it crashed.

Library(mclust)
Dat - c(rnorm(20, mean=0, sd=0.2), rnorm(30, mean=1, sd=0.2))
Ind - Mclust(dat, 1, 5)$classification
cbind(Dat, Ind)

The error message was:

/usr/local/R-2.0.1_64bit/lib/R/bin/BATCH: line 55: 18097 Done
( echo invisible(options(echo = TRUE)); cat ${in}; echo proc.time() )
 18099 Segmentation fault  | ${R_HOME}/bin/R ${opts} ${out} 21

Can anybody help me with this?
Thanks in advance,

Tae-Hoon Chung

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


Re: [R] .gct file

2005-07-19 Thread Spencer Graves
  Try ?read.table or args(read.table).  Might skip=2 do what you want?

  spencer graves
p.s.  I routinely readLines(File, n=11) to see how many headers there 
are AND identify the sep character.  Then I 
quantile(count.fields(File, ...)) to see if all records have the same 
number of fields.  Then I call something like read.table with the 
appropriate arguments.  Then I print the first 2 or so rows of the 
result to make sure I read the file correctly.

mark salsburg wrote:

 ok so the gct file looks like this:
 
 #1.2  (version number)
 7283 19   (matrix size)
 Name Description Values
   ...  ..
 
 How can I tell R to disregard the first two lines and start reading
 the 3rd line in this gct file. I would just delete them, but I do not
 know how to open a gct. file
 
 thank you
 
 On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
 
On 7/19/2005 12:10 PM, mark salsburg wrote:

I have two files to compare, one is a regular txt file that I can read
in no prob.

The other is a .gct file (How do I read in this one?)

I tried a simple

read.table(data.gct, header = T)

How do you suggest reading in this file??


.gct is not a standard filename extension.  You need to know what is in
that file.  Where did you get it?  What program created it?

Chances are the easiest thing to do is to get the program that created
it to export in a well known format, e.g. .csv.

Duncan Murdoch

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

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

[EMAIL PROTECTED]
www.pdf.com http://www.pdf.com
Tel:  408-938-4420
Fax: 408-280-7915

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


[R] Code Verification

2005-07-19 Thread pantd
Hi R Users
I have a code which I am running for my thesis work. Just want to make sure that
its ok. Its a t test I am conducting between two gamma distributions with
different shape parameters.

the code looks like:

sink(a1.txt);

for (i in 1:1000)
{
x-rgamma(40, 2.5, 10)  # n = 40, shape = 2.5, Scale = 10
y-rgamma(40, 2.8, 10)  # n = 40, shape = 2.8, Scale = 10
z-t.test(x, y)
print(z)
}


I will appreciate it if someone could tell me if its alrite or not.

thanks

-dev

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


[R] integrate fails with errors

2005-07-19 Thread [EMAIL PROTECTED]
Hi all,
i'm new to R,
I need to modelize in R a statistic algorithm,
This algo use Weibull, normal law, linear regression, normalisation, root mean 
square, to find eta and beta fitting the weibull model (to analyse few results) 
and further when we will get more information apply bayes model .

the problem is when When i try to integrate it fails with errors.

by the way i like to integrate something like that :

beta0,eta0,n are initialized as single integer, temp is a 1 dimension array 
containing 9 integer.

integrate(function(beta) 
((beta/(eta0)^beta)^n)*prod(temp^(1-beta)*exp(-sum(temp^beta)/(eta^beta)))*(1/(sqrt(2*pi))*exp(((beta-beta0)^2,0,Inf)

R says : The longuest object isn't a multiple of the shortest object in 
temp^beta , the same for temp^(1-beta).

I don't understand why R fails to take calculate each temp^beta then sum it 
over again for each beta values.

don't know if my post is explicit my head is burned out with these things today

thks.


// Webmail Oreka : http://www.oreka.com


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


Re: [R] .gct file

2005-07-19 Thread Marc Schwartz (via MN)
On Tue, 2005-07-19 at 13:16 -0400, mark salsburg wrote:
 ok so the gct file looks like this:
 
 #1.2  (version number)
 7283 19   (matrix size)
 Name Description Values
   ...  ..
 
 How can I tell R to disregard the first two lines and start reading
 the 3rd line in this gct file. I would just delete them, but I do not
 know how to open a gct. file
 
 thank you
 
 On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
  On 7/19/2005 12:10 PM, mark salsburg wrote:
   I have two files to compare, one is a regular txt file that I can read
   in no prob.
  
   The other is a .gct file (How do I read in this one?)
  
   I tried a simple
  
   read.table(data.gct, header = T)
  
   How do you suggest reading in this file??
  
  
  .gct is not a standard filename extension.  You need to know what is in
  that file.  Where did you get it?  What program created it?
  
  Chances are the easiest thing to do is to get the program that created
  it to export in a well known format, e.g. .csv.
  
  Duncan Murdoch


The above would be consistent with the info in my reply.

I guess if the format is consistent, as per Mark's example above, you
can use:

read.table(data.gct, skip = 2, header = TRUE)

which will start by skipping the first two lines and then reading in the
header row and then the data.

See ?read.table

HTH,

Marc Schwartz

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


Re: [R] .gct file

2005-07-19 Thread mark salsburg
This is all extremely helpful.

The data turns out is a little atypical, the columns are tab-delemited
except for the description columns


DATA1.gct looks like this

#1.2
23 3423
NAME DESCRIPTION VALUE
gene1 a protein inducer 1123
.  . ..

How do I get R to read the data as tab delemited, but read in the 2nd
coloumn as one value based on the quotation marks..

thanks..

On 7/19/05, Marc Schwartz (via MN) [EMAIL PROTECTED] wrote:
 On Tue, 2005-07-19 at 13:16 -0400, mark salsburg wrote:
  ok so the gct file looks like this:
 
  #1.2  (version number)
  7283 19   (matrix size)
  Name Description Values
    ...  ..
 
  How can I tell R to disregard the first two lines and start reading
  the 3rd line in this gct file. I would just delete them, but I do not
  know how to open a gct. file
 
  thank you
 
  On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
   On 7/19/2005 12:10 PM, mark salsburg wrote:
I have two files to compare, one is a regular txt file that I can read
in no prob.
   
The other is a .gct file (How do I read in this one?)
   
I tried a simple
   
read.table(data.gct, header = T)
   
How do you suggest reading in this file??
   
  
   .gct is not a standard filename extension.  You need to know what is in
   that file.  Where did you get it?  What program created it?
  
   Chances are the easiest thing to do is to get the program that created
   it to export in a well known format, e.g. .csv.
  
   Duncan Murdoch
 
 
 The above would be consistent with the info in my reply.
 
 I guess if the format is consistent, as per Mark's example above, you
 can use:
 
 read.table(data.gct, skip = 2, header = TRUE)
 
 which will start by skipping the first two lines and then reading in the
 header row and then the data.
 
 See ?read.table
 
 HTH,
 
 Marc Schwartz
 
 


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


Re: [R] using argument names (of indeterminate number) within a function

2005-07-19 Thread Gabor Grothendieck
On 7/19/05, Dirk Enzmann [EMAIL PROTECTED] wrote:
 Although I tried to find an answer in the manuals and archives, I cannot
 solve this (please excuse that my English and/or R programming skills
 are not good enough to state my problem more clearly):
 
 I want to write a function with an indeterminate (not pre-defined)
 number of arguments and think that I should use the ... construct and
 the match.call() function. The goal is to write a function that (among
 other things) uses cbind() to combine a not pre-defined number of
 vectors specified in the function call. For example, if my vectors are
 x1, x2, x3, ... xn, within the function I want to use cbind(x1, x2) or
 cbind(x1, x3, x5) or ... depending on the vector names I use in the
 funcion call. Additionally, the function has other arguments.
 
 In the archives I found the following thread (followed by Marc Schwartz)
 
 http://finzi.psych.upenn.edu/R/Rhelp02a/archive/15186.html
 [R] returning argument names from Peter Dalgaard BSA on 2003-04-10 (stdin)
 
 that seems to contain the solution to my problem, but I am stuck because
  sapply(match.call()[-1], deparse) gives me a vector of strings and I
 don't know how to use the names in this vector in the cbind() function.
 
 Up to now my (clearly deficit) function looks like:
 
 test - function(..., mvalid=1)
 {
   args = sapply(match.call()[-1], deparse)
 # and here, I don't know how the vector names in args
 # can be used in the cbind() function to follow:
 #
 # temp - cbind( ???
   if (mvalid  1)
   {
 #  here it goes on
   }
 }
 
 Ultimately, I want that the function can be called like
 test(x1,x2,mvalid=1)
 or
 test(x1,x3,x5,mavlid=2)
 and that within the function
 cbind(x1,x2)
 or cbind(x1,x3,x5)
 will be used.
 

If you just need to pass them to cbind then just use cbind(...), e.g.

  test - function(..., m) if (m  1) cbind(...) else m

otherwise, use list(...) as shown by a previous answer or here

  test2 - function(..., m) if (m  1) length(list(...)) else m

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


Re: [R] .gct file

2005-07-19 Thread Marc Schwartz (via MN)
For the TAB delimited columns, adjust the 'sep' argument to:

read.table(data.gct, skip = 2, header = TRUE, sep = \t)

The 'quote' argument is by default:

quote = \'

which should take care of the quoted strings and bring them in as a
single value.

The above presumes that the header row is also TAB delimited. If not,
you may have to set 'skip = 3' to skip over the header row and manually
set the column names.

HTH,

Marc Schwartz


On Tue, 2005-07-19 at 13:52 -0400, mark salsburg wrote:
 This is all extremely helpful.
 
 The data turns out is a little atypical, the columns are tab-delemited
 except for the description columns
 
 
 DATA1.gct looks like this
 
 #1.2
 23 3423
 NAME DESCRIPTION VALUE
 gene1 a protein inducer 1123
 .  . ..
 
 How do I get R to read the data as tab delemited, but read in the 2nd
 coloumn as one value based on the quotation marks..
 
 thanks..
 
 On 7/19/05, Marc Schwartz (via MN) [EMAIL PROTECTED] wrote:
  On Tue, 2005-07-19 at 13:16 -0400, mark salsburg wrote:
   ok so the gct file looks like this:
  
   #1.2  (version number)
   7283 19   (matrix size)
   Name Description Values
     ...  ..
  
   How can I tell R to disregard the first two lines and start reading
   the 3rd line in this gct file. I would just delete them, but I do not
   know how to open a gct. file
  
   thank you
  
   On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
On 7/19/2005 12:10 PM, mark salsburg wrote:
 I have two files to compare, one is a regular txt file that I can read
 in no prob.

 The other is a .gct file (How do I read in this one?)

 I tried a simple

 read.table(data.gct, header = T)

 How do you suggest reading in this file??

   
.gct is not a standard filename extension.  You need to know what is in
that file.  Where did you get it?  What program created it?
   
Chances are the easiest thing to do is to get the program that created
it to export in a well known format, e.g. .csv.
   
Duncan Murdoch
  
  
  The above would be consistent with the info in my reply.
  
  I guess if the format is consistent, as per Mark's example above, you
  can use:
  
  read.table(data.gct, skip = 2, header = TRUE)
  
  which will start by skipping the first two lines and then reading in the
  header row and then the data.
  
  See ?read.table
  
  HTH,
  
  Marc Schwartz
  
  
 

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


Re: [R] Michaelis-menten equation

2005-07-19 Thread Peter Dalgaard
Peter Dalgaard [EMAIL PROTECTED] writes:

 Chun-Ying Lee [EMAIL PROTECTED] writes:
 
  Dear R users:
 I encountered difficulties in michaelis-menten equation. I found 
  that when I use right model definiens, I got wrong Km vlaue, 
  and I got right Km value when i use wrong model definiens. 
  The value of Vd and Vmax are correct in these two models. 
 
 How do you know what the correct value is? Are you sure that the other
 values are right?
 
 I'm a bit rusty on MM, but are you sure your right model is right?
 Try doing a dimensional analysis on the ODE. I kind of suspect that
 Vd is entering in the wrong way. Since you're dealing in
 concentrations, should it enter at all (except via the conc. at time
 0, of course)?
 
 Not knowing the context, I can't be quite sure, but generally, I'd
 expect Vm*Km/(Km+y) to be the reaction rate, so that Vm is the maximum
 rate, attained when y is zero and Km is the conc. at half-maximum
 rate. This doesn't look quit like what you have. 

Hmm, sorry, no. I'm talking through a hole in my head there.

Vm*y/(Km+y) makes OK sense. Vm is what you get for large y - passing
from 1st order to 0th order kinetics. However, looking at the data

 plot(PKindex)
 abline(lm(conc~time,data=PKindex))

shows that they are pretty much on a straight line, i.e. you are 
in the domain of 0-order kinetics. So why are you expecting the rate
of decrease to have changed by roughly 3/4 (from 2/3*Vm/Vd at y=2*Km
to 1/2*Vm/Vd at y=Km when you reach 4.67)??
  
  #-right model definiens
  PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
 conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
  mm.model - function(time, y, parms) { 
 dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]+y[1]) 
 list(dCpdt)}
  Dose-300
  modfun - function(time,Vm,Km,Vd) { 
 out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
rtol=1e-8,atol=1e-8)
out[,2] } 
  objfun - function(par) { 
 out - modfun(PKindex$time,par[1],par[2],par[3]) 
 sum((PKindex$conc-out)^2) } 
  fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
  print(fit$par)
  [1] 10.0390733  0.1341544 34.9891829  #--Km=0.1341544,wrong value--
  
  
  #-wrong model definiens
  #-Km should not divided by Vd--
  PKindex-data.frame(time=c(0,1,2,4,6,8,10,12,16,20,24),
 conc=c(8.57,8.30,8.01,7.44,6.88,6.32,5.76,5.20,4.08,2.98,1.89))
  mm.model - function(time, y, parms) { 
 dCpdt - -(parms[Vm]/parms[Vd])*y[1]/(parms[Km]/parms[Vd]+y[1]) 
 list(dCpdt)}
  Dose-300
  modfun - function(time,Vm,Km,Vd) { 
  out - lsoda(Dose/Vd,time,mm.model,parms=c(Vm=Vm,Km=Km,Vd=Vd),
  rtol=1e-8,atol=1e-8)
 out[,2] 
  } 
  objfun - function(par) { 
  out - modfun(PKindex$time,par[1],par[2],par[3]) 
  sum((PKindex$conc-out)^2)} 
  fit - optim(c(10,1,80),objfun, method=Nelder-Mead)
  print(fit$par)
  [1] 10.038821  4.690267 34.989239  #--Km=4.690267,right value--
  
  What did I do wrong, and how to fix it?
  Any suggestions would be greatly appreciated.
  Thanks in advance!!
  
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide! 
  http://www.R-project.org/posting-guide.html
 
 -- 
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
 ~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 

-- 
   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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

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


Re: [R] .gct file

2005-07-19 Thread Marc Schwartz
On Tue, 2005-07-19 at 13:08 -0500, Marc Schwartz (via MN) wrote:
 For the TAB delimited columns, adjust the 'sep' argument to:
 
 read.table(data.gct, skip = 2, header = TRUE, sep = \t)
 
 The 'quote' argument is by default:
 
 quote = \'
 
 which should take care of the quoted strings and bring them in as a
 single value.
 
 The above presumes that the header row is also TAB delimited. If not,
 you may have to set 'skip = 3' to skip over the header row and manually
 set the column names.

One correction. If the final para applies and you need to use 'skip =
3', you would also need to leave out the 'header = TRUE' argument, which
defaults to FALSE.

Marc

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


Re: [R] .gct file

2005-07-19 Thread Marc Schwartz (via MN)
On Tue, 2005-07-19 at 13:08 -0500, Marc Schwartz (via MN) wrote:
 For the TAB delimited columns, adjust the 'sep' argument to:
 
 read.table(data.gct, skip = 2, header = TRUE, sep = \t)
 
 The 'quote' argument is by default:
 
 quote = \'
 
 which should take care of the quoted strings and bring them in as a
 single value.
 
 The above presumes that the header row is also TAB delimited. If not,
 you may have to set 'skip = 3' to skip over the header row and manually
 set the column names.

One correction. If the final para applies and you need to use 'skip =
3', you would also need to leave out the 'header = TRUE' argument, which
defaults to FALSE.

Marc

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


[R] CPU Usage with R 2.1.0 in Windows

2005-07-19 Thread Greene, Michael

Hi,

I'm using a fairly simple HP Compaq desktop PC running Windows 2K.  When
running a large process in R, the process RGUI.exe will never exceed 50%
of the CPU usage.

The program used to be able to use more of the computer, but does not now.
I don't believe this is a multiple processor machine.

Can anyone give any advice on how to solve the problem?  

Thanks,

Michael Greene

Product Management
Plymouth Rock Assurance Corp
617-951-1682

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


Re: [R] ROC curve with survival data

2005-07-19 Thread Frank E Harrell Jr
SUBIRANA CACHINERO, ISAAC wrote:
 Hi everyone, 
 
 I am doing 5 years mortality predictive index score with survival analysis 
 using a Cox proportional hazard model where I have a continous predictive 
 variable and a right censored response which is the mortality, and the 
 individuals were followed a maximum of 7 years.
 
 I'd like to asses the discrimination ability of survival analysis Cox model 
 by computing a ROC curve and area under the curve for a fixed time (5 years), 
 taking into account that the response is not binary but right censored.
 
 So, is there a function that computes a ROC curve under right censored 
 survival data?
 
 Thank you in advance.

I don't find ROC curves themselves very useful but the area under them 
is useful and is a simple translation of the Somers' Dxy rank 
correlation between predicted survival probability (or anything 
monotonically related to it just is log hazard) and observed survival 
time.  Dxy = 2*(C-.5) where C is the concordance index, a generalization 
of ROC area not requiring choosing a specific time point.  You can get 
this in the rcorr.cens function in the Hmisc package.

Frank

 
  
 
  
 
 Isaac Subirana: [EMAIL PROTECTED]
 
 Institut Municipal d'Investigació Mèdica (IMIM),
 
 Barcelona (Spain)
 
  
 
 
   [[alternative HTML version deleted]]
 
 
 
 
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


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

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


[R] Taking the derivative of a quadratic B-spline

2005-07-19 Thread James McDermott
Hello,

I have been trying to take the derivative of a quadratic B-spline
obtained by using the COBS library.  What I would like to do is
similar to what one can do by using

fit-smooth.spline(cdf)
xx-seq(-10,10,.1)
predict(fit, xx, deriv = 1)

The goal is to fit the spline to data that is approximating a
cumulative distribution function (e.g. in my example, cdf is a
2-column matrix with x values in column 1 and the estimate of the cdf
evaluated at x in column 2) and then take the first derivative over a
range of values to get density estimates.

The reason I don't want to use smooth.spline is that there is no way
to impose constraints (e.g. =0, =1, and monotonicity) as there is
with COBS.  However, since COBS doesn't have the 'deriv =' option, the
only way I can think of doing it with COBS is to evaluate the
derivatives numerically.

Regards,
Jim McDermott

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


Re: [R] CPU Usage with R 2.1.0 in Windows

2005-07-19 Thread Doran, Harold
Dear Michael:

Why is it a problem that R is not using more CPU space than it seems to
need? 

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Greene, Michael
Sent: Tuesday, July 19, 2005 2:29 PM
To: '[EMAIL PROTECTED]'
Subject: [R] CPU Usage with R 2.1.0 in Windows


Hi,

I'm using a fairly simple HP Compaq desktop PC running Windows 2K.  When
running a large process in R, the process RGUI.exe will never exceed
50% of the CPU usage.

The program used to be able to use more of the computer, but does not
now.
I don't believe this is a multiple processor machine.

Can anyone give any advice on how to solve the problem?  

Thanks,

Michael Greene

Product Management
Plymouth Rock Assurance Corp
617-951-1682

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

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


Re: [R] Taking the derivative of a quadratic B-spline

2005-07-19 Thread Duncan Murdoch
On 7/19/2005 2:53 PM, James McDermott wrote:
 Hello,
 
 I have been trying to take the derivative of a quadratic B-spline
 obtained by using the COBS library.  What I would like to do is
 similar to what one can do by using
 
 fit-smooth.spline(cdf)
 xx-seq(-10,10,.1)
 predict(fit, xx, deriv = 1)
 
 The goal is to fit the spline to data that is approximating a
 cumulative distribution function (e.g. in my example, cdf is a
 2-column matrix with x values in column 1 and the estimate of the cdf
 evaluated at x in column 2) and then take the first derivative over a
 range of values to get density estimates.
 
 The reason I don't want to use smooth.spline is that there is no way
 to impose constraints (e.g. =0, =1, and monotonicity) as there is
 with COBS.  However, since COBS doesn't have the 'deriv =' option, the
 only way I can think of doing it with COBS is to evaluate the
 derivatives numerically.

Numerical estimates of the derivatives of a quadratic should be easy to 
obtain accurately.  For example, if the quadratic ax^2 + bx + c is 
defined on [-1, 1], then the derivative 2ax + b, has 2a = f(1) - f(0) + 
f(-1), and b = (f(1) - f(-1))/2.

You should be able to generalize this to the case where the spline is 
quadratic between knots k1 and k2 pretty easily.

Duncan Murdoch

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


[R] Question about creating unique factor labels with the factor function

2005-07-19 Thread Gregory Gentlemen
Hi guys,
 
I ran into a problem of not being able to create unique labels when creating a 
factor. Consider an example below:
 
 hb - factor(c(1,1,1,2,2,2,3,3,3), levels=c(1,2,3),labels=c(1,1,2))
 hb
[1] 1 1 1 1 1 1 2 2 2
Levels: 1 1 2
 unique(hb)
[1] 1 1 2
Levels: 1 1 2
 
How come there are three unique levels, I thought this would only create one 
unique level?

 unique(as.ordered(hb))
[1] 1 2
Levels: 1  1  2
 
Is as.ordered the only solution?
 
Thanks in advance,
Greg
 

__



[[alternative HTML version deleted]]

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


Re: [R] CPU Usage with R 2.1.0 in Windows

2005-07-19 Thread Lukasz Komsta
Dnia 2005-07-19 20:28, Użytkownik Greene, Michael napisał:

 Hi,
 
 I'm using a fairly simple HP Compaq desktop PC running Windows 2K.  When
 running a large process in R, the process RGUI.exe will never exceed 50%
 of the CPU usage.

If you have hyperthreading, R catches only one virtual processor (from
two available), being not able to exceed half of total power (100% of
one only). If you want to use full power, you should turn hyperthreading
off, if your BIOS supports such option.

Regards,

-- 
Lukasz Komsta
Department of Medicinal Chemistry
Medical University of Lublin
6 Chodzki, 20-093 Lublin, Poland
Fax +48 81 7425165

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

Re: [R] Taking the derivative of a quadratic B-spline

2005-07-19 Thread James McDermott
I wish it were that simple (perhaps it is and I am just not seeing
it).  The output from cobs( ) includes the B-spline coefficients and
the knots.  These coefficients are not the same as the a, b, and c
coefficients in a quadratic polynomial.  Rather, they are the
coefficients of the quadratic B-spline representation of the fitted
curve.  I need to evaluate a linear combination of basis functions and
it is not clear to me how to accomplish this easily.  I was hoping to
find an alternative way of getting the derivatives.

Jim McDermott

On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
 On 7/19/2005 2:53 PM, James McDermott wrote:
  Hello,
 
  I have been trying to take the derivative of a quadratic B-spline
  obtained by using the COBS library.  What I would like to do is
  similar to what one can do by using
 
  fit-smooth.spline(cdf)
  xx-seq(-10,10,.1)
  predict(fit, xx, deriv = 1)
 
  The goal is to fit the spline to data that is approximating a
  cumulative distribution function (e.g. in my example, cdf is a
  2-column matrix with x values in column 1 and the estimate of the cdf
  evaluated at x in column 2) and then take the first derivative over a
  range of values to get density estimates.
 
  The reason I don't want to use smooth.spline is that there is no way
  to impose constraints (e.g. =0, =1, and monotonicity) as there is
  with COBS.  However, since COBS doesn't have the 'deriv =' option, the
  only way I can think of doing it with COBS is to evaluate the
  derivatives numerically.
 
 Numerical estimates of the derivatives of a quadratic should be easy to
 obtain accurately.  For example, if the quadratic ax^2 + bx + c is
 defined on [-1, 1], then the derivative 2ax + b, has 2a = f(1) - f(0) +
 f(-1), and b = (f(1) - f(-1))/2.
 
 You should be able to generalize this to the case where the spline is
 quadratic between knots k1 and k2 pretty easily.
 
 Duncan Murdoch


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


Re: [R] Taking the derivative of a quadratic B-spline

2005-07-19 Thread Huntsinger, Reid
The derivative of a quadratic B-spline is the centered finite difference of
a linear B-spline, so if you have access to the underlying coefficients of
the B-spline expansion you can do this easily. I believe the coefficients
are passed as the $coef component of the return value.

Reid Huntsinger

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of James McDermott
Sent: Tuesday, July 19, 2005 2:54 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Taking the derivative of a quadratic B-spline


Hello,

I have been trying to take the derivative of a quadratic B-spline
obtained by using the COBS library.  What I would like to do is
similar to what one can do by using

fit-smooth.spline(cdf)
xx-seq(-10,10,.1)
predict(fit, xx, deriv = 1)

The goal is to fit the spline to data that is approximating a
cumulative distribution function (e.g. in my example, cdf is a
2-column matrix with x values in column 1 and the estimate of the cdf
evaluated at x in column 2) and then take the first derivative over a
range of values to get density estimates.

The reason I don't want to use smooth.spline is that there is no way
to impose constraints (e.g. =0, =1, and monotonicity) as there is
with COBS.  However, since COBS doesn't have the 'deriv =' option, the
only way I can think of doing it with COBS is to evaluate the
derivatives numerically.

Regards,
Jim McDermott

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

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


Re: [R] Taking the derivative of a quadratic B-spline

2005-07-19 Thread Duncan Murdoch
On 7/19/2005 3:34 PM, James McDermott wrote:
 I wish it were that simple (perhaps it is and I am just not seeing
 it).  The output from cobs( ) includes the B-spline coefficients and
 the knots.  These coefficients are not the same as the a, b, and c
 coefficients in a quadratic polynomial.  Rather, they are the
 coefficients of the quadratic B-spline representation of the fitted
 curve.  I need to evaluate a linear combination of basis functions and
 it is not clear to me how to accomplish this easily.  I was hoping to
 find an alternative way of getting the derivatives.

I don't know COBS, but doesn't predict just evaluate the B-spline?  The 
point of what I posted is that the particular basis doesn't matter if 
you can evaluate the quadratic at 3 points.

Duncan Murdoch

 
 Jim McDermott
 
 On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
 On 7/19/2005 2:53 PM, James McDermott wrote:
  Hello,
 
  I have been trying to take the derivative of a quadratic B-spline
  obtained by using the COBS library.  What I would like to do is
  similar to what one can do by using
 
  fit-smooth.spline(cdf)
  xx-seq(-10,10,.1)
  predict(fit, xx, deriv = 1)
 
  The goal is to fit the spline to data that is approximating a
  cumulative distribution function (e.g. in my example, cdf is a
  2-column matrix with x values in column 1 and the estimate of the cdf
  evaluated at x in column 2) and then take the first derivative over a
  range of values to get density estimates.
 
  The reason I don't want to use smooth.spline is that there is no way
  to impose constraints (e.g. =0, =1, and monotonicity) as there is
  with COBS.  However, since COBS doesn't have the 'deriv =' option, the
  only way I can think of doing it with COBS is to evaluate the
  derivatives numerically.
 
 Numerical estimates of the derivatives of a quadratic should be easy to
 obtain accurately.  For example, if the quadratic ax^2 + bx + c is
 defined on [-1, 1], then the derivative 2ax + b, has 2a = f(1) - f(0) +
 f(-1), and b = (f(1) - f(-1))/2.
 
 You should be able to generalize this to the case where the spline is
 quadratic between knots k1 and k2 pretty easily.
 
 Duncan Murdoch


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


Re: [R] extracting row means from a list

2005-07-19 Thread Andy Bunn
I think about half of my question in R can be solved with a judicious
do.call.

Thanks, Andy

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


[R] Using BRugs, in FUN: .C(..) 'type' not real

2005-07-19 Thread Seth Pruitt
To All,

I am using the BRugs package. In running the meta file BRugsFit() with a 
syntactically correct model .txt file, I see the message:

Error in FUN(X[[1]], ...) : .C(..): 'type' must be real for this format

I haven't found information on this kind of error either here or on the 
OpenBUGS/WinBUGS/BRugs bulletin boards. If there is a better place to post 
this message, please let me know.

I am running Windows XP Professional (build 2600) Service Pack 2.0 on Intel 
Centrino, and using R version 2.1.1.

Thank you,

-- 
Seth Pruitt
Department of Economics
University of California, San Diego
[EMAIL PROTECTED]
http://dss.ucsd.edu/~sjpruitt

[[alternative HTML version deleted]]

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


[R] A warning message when using mix package.

2005-07-19 Thread Li, Jia
Dear all, 

I just start using package mix in R and am inexperienced with it. When I
ran my program several times, I sometimes got the warning message, and
sometimes not : 

Warning message:

Loglik converged before variable 2 ; beta may be infinite. in: fitter(X,
Y, strats, offset, init, control, weights = weights, 

I cannot find what is wrong with my program. Does anyone have an idea
with it and whether this would affect my results? Any help would be
appreciated.

Thanks!

By the way, I copy a part of my program and hope it would help,

beta - c(0, 0,0,0,0) # 'True' coefficients

#make a matrix of covariates---#

Y=cbind(z4,z5,z1,z2,z3) 

s-prelim.mix(Y,2)

#run MI--#

MI-vector(list,10)

fit.model.mi-vector(list,10)

rngseed(1234567) # set random number generator seed

for (i in 1:10){

thetahat - em.mix(s)

newtheta - da.mix(s, thetahat, steps=2000, showits=TRUE) 

MI[[i]] - imp.mix(s, newtheta)

}


[[alternative HTML version deleted]]

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


[R] sciviews installation

2005-07-19 Thread ghjmora g mail
Hello


1. a few months ago, I had sciviews working fine with R (rw2001) under 
windows XP
2. now, upgrading to rw2011, the stuff seems fine (every package 
installed),but I find a conflict when launching sciviews:
- it runs, apparently
- but when I try to work (import/export In: text for instance), it 
asks for Rcmdr (Would you like to install it now?)

3. Rcmdr is already installed (with all dependencies) and works well 
when called directly in R gui
4. and it's impossible to make it reconized or to install it under sciviews


I have all the latest packages, and I am going to get mad.

what do you suggest to solve my problem ?

Thanks

Georges Moracchini

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


Re: [R] initial points for arms in package HI

2005-07-19 Thread plummer
Quoting Christoph Buser [EMAIL PROTECTED]:

 Dear R-users

 I have a problem choosing initial points for the function arms()
 in the package HI
 I intend to implement a Gibbs sampler and one of my conditional
 distributions is nonstandard and not logconcave.
 Therefore I'd like to use arms.

 But there seem to be a strong influence of the initial point
 y.start. To show the effect I constructed a demonstration
 example. It is reproducible without further information.

 Please note that my target density is not logconcave.

 Thanks for all comments or ideas.

 Christoph Buser

Dear Christoph,

There is a Metropolis step at each iteration of the ARMS sampler, in
which it may choose to reject the proposed move to a new point and stick
at the current point (This is what the M in ARMS stands for)  If you
do repeated calls to arms with the same starting point, then the
iterations where the Metropolis step rejects a move will create a spike
in the sample density at your initial value. If you use a uniform random
starting point, then your sample density will be a mixture of the
target distribution (Metropolis accepts move) and a uniform distribution
(Metropolis rejects move).

You should be doing something like this:

res1 - arms(runif(1,0,100), logDichteGam, function(x) (x0)(x100), 1)
for(i in 2:1000)
  res1[i] - arms(res1[i-1], logDichteGam, function(x) (x0)(x100), 1)

i.e., using each sampled point as the starting value for the next
iteration.  The sequence of values in res1 will then be a correlated
sample from the given distribution:

acf(res1)

The bottom line is that you can't use ARMS to draw a single sample
from a non-log-concave density.

If you are still worried about using ARMS, you can verify your results
using the random walk Metropolis sampler (MCMCmetrop1R) in the package
MCMCpack.

Martyn

 ## R Code:

 library(HI)
 ## parameter for the distribution
 para - 0.1

 ## logdensity
 logDichteGam - function(x, u = para, v = para) {
   -(u*x + v*1/x) - log(x)
 }
 ## density except for the constant
 propDichteGam - function(x, u = para, v = para) {
   exp(-(u*x + v*1/x) - log(x))
 }
 ## calculating the constant
 (c - integrate(propDichteGam, 0, 1000, rel.tol = 10^(-12))$value)
 ## density
 DichteGam - function(x, u = para, v = para) {
   exp(-(u*x + v*1/x) - log(x))/c
 }

 ## calculating 1000 values by repeating a call of arms (this would
 ## be the situation in an Gibbs Sample. Of course in a Gibbs sampler
 ## the distribution would change. This is only for demonstration
 res1 - NULL
 for(i in 1:1000)
   res1[i] - arms(runif(1,0,100), logDichteGam, function(x) (x0)(x100), 1)

 ## Generating a sample of thousand observations with 1 call of arms
 res2 - arms(runif(1,0,100), logDichteGam, function(x) (x0)(x100), 1000)

 ## Plot of the samples
 mult.fig(4)
 plot(res1, log = y)
 plot(res2, log = y)
 hist(res1, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
  ylim = c(0,1))
 curve(DichteGam, 0,4, add = TRUE, col = 2)
 hist(res2, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
  ylim = c(0,1))
 curve(DichteGam, 0,4, add = TRUE, col = 2)


 ## If we repeat the procedure, using the fix intial value 1,
 ## the situation is even worse
 res3 - NULL
 for(i in 1:1000)
   res3[i] - arms(1, logDichteGam, function(x) (x0)(x100), 1)

 ## Generating a sample of thousand observations with 1 call of arms
 res4 - arms(1, logDichteGam, function(x) (x0)(x100), 1000)

 ## Plot of the samples
 par(mfrow = c(2,2))
 plot(res3, log = y)
 plot(res4, log = y)
 hist(res3, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
  ylim = c(0,1))
 curve(DichteGam, 0,4, add = TRUE, col = 2)
 hist(res4, freq = FALSE, xlim = c(0,4), breaks = seq(0,100,by = 0.1),
  ylim = c(0,1))
 curve(DichteGam, 0,4, add = TRUE, col = 2)


 ## If I generate the sample in a for-loop (one by one) I do not
 ## get the correct density. But this is exactly the situation in
 ## my Gibbs Sampler. Therfore I am concerned about the correct
 ## application of arms




---
This message and its attachments are strictly confidential. ...{{dropped}}

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


[R] Is it possible to create highly customized report in *.xls format by using R/S+?

2005-07-19 Thread Wensui Liu
I remember in one slide of Prof. Ripley's presentation overhead, he
said the most popular data analysis software is excel.

So is there any resource or tutorial on this topic? 

Thank you so much!

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


[R] Problems with date-format (R 2.1.1 + chron)

2005-07-19 Thread Carsten Steinhoff
Hello,
 
today I've updated on the newest R-Version. But sadly a function I needed
didnt want to work:
The input is e.g.
 
days(as.Date(21-07-2005,%d-%m-%y))

the error is: Fehler in Math.Date(dts): floor nicht definiert für Date
Objekte
(Error in Math.Date(dts): floor not defined for date objects)

Same for year. Only months gives me the correct output.
In Version 2.01 it worked very well, with the same chron library.
Whats wrong ?
 
Carsten

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


[R] Indexing within scan

2005-07-19 Thread Justin Rhodes
Dear R-help subscribers,

Can some one please help me figure out how to write code that will allow me to 
use a for loop to scan a number of files one by one, and then save a summary of 
each file as the for loop progresses.  For example I have 24 files named a1 
through a24, and I want to do something like:


results-numeric(24)
for (i in 1:24)
{
 
p-scan(ai.txt)  # where the i is an index for each of the 24 files

results[i]-mean(p)
}

Thanks for any help,





Justin Rhodes
Behavioral Neuroscience
Oregon Health  Science University
VA Medical Center (R  D 12)
3710 SW US Veterans Hospital Rd
Portland, OR  97239
Phone: (503) 220-8262 extn 54392
Fax: (503) 721-1029
E-mail: [EMAIL PROTECTED]

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


Re: [R] Indexing within scan

2005-07-19 Thread Peter Dalgaard
Justin Rhodes [EMAIL PROTECTED] writes:

 Dear R-help subscribers,
 
 Can some one please help me figure out how to write code that will allow me 
 to use a for loop to scan a number of files one by one, and then save a 
 summary of each file as the for loop progresses.  For example I have 24 files 
 named a1 through a24, and I want to do something like:
 
 
 results-numeric(24)
 for (i in 1:24)
 {
  
 p-scan(ai.txt)  # where the i is an index for each of the 24 files
 
 results[i]-mean(p)
 }
 
 Thanks for any help,
 

paste(a, i, .txt, sep=) 

should get you in the right direction. However, I don't think you
*really* want a for loop. Use vectorization and sapply instead:

names - paste(a, 1:24, .txt, sep=)
results - sapply(names, function(n) mean(scan(n)) )

-- 
   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
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

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


Re: [R] extracting row means from a list

2005-07-19 Thread Liaw, Andy
Justone more comment in addition to Sundar's solution:

If these are all numeric matrices, I would read them into R
as such, instead of data frames.  Actually, I would read them
all into a 3-dimensional array (2000 x 6 x # of files).
Assuming you have such an array, then you can do something
like:

 ## get your example into an array: need the abind package.
 bar - do.call(abind, c(lapply(foo, as.matrix), along=3))
 str(bar)
 num [1:100, 1:6, 1:10] -0.78981  0.31939 -0.00819  1.59346  1.20498 ...
 - attr(*, dimnames)=List of 3
  ..$ : chr [1:100] 1 2 3 4 ...
  ..$ : chr [1:6] x1 x2 x3 x4 ...
  ..$ : NULL
 str(m - rowMeans(bar, dims=2))
 num [1:100, 1:6] -0.4401  0.5463 -0.0572 -0.1314  0.5177 ...
 - attr(*, dimnames)=List of 2
  ..$ : chr [1:100] 1 2 3 4 ...
  ..$ : chr [1:6] x1 x2 x3 x4 ...

Andy


 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of Andy Bunn
 Sent: Tuesday, July 19, 2005 3:59 PM
 To: Sundar Dorai-Raj
 Cc: R-Help
 Subject: Re: [R] extracting row means from a list
 
 
 I think about half of my question in R can be solved with a judicious
 do.call.
 
 Thanks, Andy
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! 
 http://www.R-project.org/posting-guide.html
 
 


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


[R] deriv - accessing numeric output for gradient

2005-07-19 Thread Rotella, Jay
Hi,
I am interested in using the numeric output from the gradient attribute of 
deriv's output in subsequent analyses.
But, I have so far been unable to determine how to do so.

I will use the example from the deriv help to illustrate.

 ## function with defaulted arguments:
(fx - deriv(y ~ b0 + b1 * 2^(-x/th), c(b0, b1, th),
  function(b0, b1, th, x = 1:7){} ) )
 fx(2,3,4)

This yields
[1] 4.522689 4.121320 3.783811 3.50 3.261345 3.060660 2.891905
attr(,gradient)
 b0b1th
[1,]  1 0.8408964 0.1092872
[2,]  1 0.7071068 0.1837984
[3,]  1 0.5946036 0.2318331
[4,]  1 0.500 0.2599302
[5,]  1 0.4204482 0.2732180
[6,]  1 0.3535534 0.2756976
[7,]  1 0.2973018 0.2704720

I would greatly appreciate it if anyone could tell me how to convert the 
numbers listed under b0, b1, and th into a matrix.

Thanks!

Jay Rotella
Ecology Department
Montana State University
Bozeman, MT 59717







[[alternative HTML version deleted]]

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


[R] deriv - accessing numeric output listed under gradient attribute

2005-07-19 Thread Jay Rotella
Hi,
I am interested in using the numeric output from the gradient attribute of 
deriv's output in subsequent analyses.
But, I have so far been unable to determine how to do so.

I will use the example from the deriv help to illustrate.

 ## function with defaulted arguments:
(fx - deriv(y ~ b0 + b1 * 2^(-x/th), c(b0, b1, th),
  function(b0, b1, th, x = 1:7){} ) )
 fx(2,3,4)

This yields
[1] 4.522689 4.121320 3.783811 3.50 3.261345 3.060660 2.891905
attr(,gradient)
 b0b1th
[1,]  1 0.8408964 0.1092872
[2,]  1 0.7071068 0.1837984
[3,]  1 0.5946036 0.2318331
[4,]  1 0.500 0.2599302
[5,]  1 0.4204482 0.2732180
[6,]  1 0.3535534 0.2756976
[7,]  1 0.2973018 0.2704720

I would greatly appreciate it if anyone could tell me how to convert the 
numbers listed under b0, b1, and th into a matrix.

Thanks!

Jay Rotella
Ecology Department
Montana State University
Bozeman, MT 59717

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



[R] Regression lines for differently-sized groups on the same plot

2005-07-19 Thread Laura M Marx
Hi there,
  I've looked through the very helpful advice about adding fitted lines to 
plots in the r-help archive, and can't find a post where someone has offered 
a solution for my specific problem.  I need to plot logistic regression fits 
from three differently-sized data subsets on a plot of the entire dataset.  
A description and code are below:
  I have an unbalanced dataset consisting of three different species (hem, 
yb, and sm), with unequal numbers of wood pieces in each species group.  I 
am trying to generate a plot that will show the size of the wood piece on 
the X axis, the probability of it having tree seedlings growing on it on the 
Y (a binomial yes or no variable), and three fitted curves showing how the 
probability of having tree seedlings changes with increasing wood piece size 
for each species.
  I have no problem generating fits using GLM, and no problem creating the 
plot.  However, if I try to add a fitted curve based only on the hem data 
subset to a plot that shows the entire dataset, I get an error message that 
the lengths of those data sets differ. Error in xy.coords(x,y) : x and y 
lengths differ.  I could see R's point -- you can't plot a regression line 
of babies born as a function of stork abundance on a graph of cherries 
produced (Y) versus rainfall (X), which for all the program knows, I'm 
trying to do.  As a temporary fix, I added NAs to the end of the hem, yb, 
and sm subsets to make them the same length as the entire dataset.  I can 
now add my fitted curves to the plot, but the lines are not connected.  That 
is, if the hem group only contains wood pieces that are 1, 4, and 10 meters 
long, the plot has an X axis that ranges from 1 to 10, but line segments for 
the hem group regression line only appear above 1, 4, and 10.  How can I fix 
this?  An ideal solution would not require me to make the hem subset of my 
data the same length as the full dataset, either (although the summaries of 
regressions with the NAs (or zeroes) added and taken away are identical).  
I'd also settle for a work-around that would have R connect the pieces of 
the curve so that I get a solid line rather than small dots and dashes where 
actual data exist.  Thanks so much for your help!
  Laura Marx
  Michigan State University, Dept. of Forestry 

#Note: hemdata has all the rows that are not hemlock species replaced with 
#NAs.
hemhem=glm(hempresence~logarea, family=binomial(logit), data=hemdata)
hemyb=glm(hempresence~logarea, family=binomial(logit), data=birchdata)
hemsm=glm(hempresence~logarea, family=binomial(logit), data=mapledata) 

attach(logreg) #logreg is the full dataset
plot(logarea, hempresence, xlab = Surface area of log (m2), 
ylab=Probability of hemlock seedling presence, type=n, font.lab=2, 
cex.lab=1.5, axes=TRUE)
lines(logarea,fitted(hemhem), lty=1, lwd=2)
lines(logarea,fitted(hemyb), lty=dashed, lwd=2)
lines(logarea,fitted(hemsm), lty=dotted, lwd=2)

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


Re: [R] deriv - accessing numeric output listed under gradient attribute

2005-07-19 Thread Sundar Dorai-Raj


Jay Rotella wrote:
 Hi,
 I am interested in using the numeric output from the gradient attribute of 
 deriv's output in subsequent analyses.
 But, I have so far been unable to determine how to do so.
 
 I will use the example from the deriv help to illustrate.
 
 
## function with defaulted arguments:
   (fx - deriv(y ~ b0 + b1 * 2^(-x/th), c(b0, b1, th),
 
   function(b0, b1, th, x = 1:7){} ) )
 
fx(2,3,4)
 
 
 This yields
 [1] 4.522689 4.121320 3.783811 3.50 3.261345 3.060660 2.891905
 attr(,gradient)
  b0b1th
 [1,]  1 0.8408964 0.1092872
 [2,]  1 0.7071068 0.1837984
 [3,]  1 0.5946036 0.2318331
 [4,]  1 0.500 0.2599302
 [5,]  1 0.4204482 0.2732180
 [6,]  1 0.3535534 0.2756976
 [7,]  1 0.2973018 0.2704720
 
 I would greatly appreciate it if anyone could tell me how to convert the 
 numbers listed under b0, b1, and th into a matrix.
 
 Thanks!
 
 Jay Rotella
 Ecology Department
 Montana State University
 Bozeman, MT 59717
 


Try:

attr(fx(2, 3, 4), gradient)

--sundar

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


Re: [R] Memory limits using read.table on Windows XP Pro

2005-07-19 Thread Latchezar Dimitrov
Hello everyone,

Would you please somebody explain me what my sin is (please see the code
and timing bellow)? And how to improve myself and the following piece of
R code? BTW, the code works.

This is R 64-bit built by myself on Sun SPARC Solaris 9 with gcc-4.0.1
(64-bit) also built by yours truly. The machine is Sun Fire-V880 with
32GB memory and 8 cpu's. Nobody else was using it.


Thank you very much

Latchezar Dimitrov

PS. Prof. Brian D. Ripley, thank you very much for not wasting your
invaluable time responding to the above. I know your answer anyway.
Sorry for making you read it though ...

  system.time(
+ haplo[i,2*j-1]-substr(as.character(geno[i,j]),1,1)
+  ,TRUE)
[1] 66.27 13.67 80.02  0.00  0.00
 ls()
[1] P geno  haplo i j lrmodel pheno
 i
[1] 1
 j
[1] 1
 dim(P)
[1]  1 125000
 str(P)
 num [1, 1:125000] 0.6188 0.0533 0.0893 0.8994 0.0316 ...
 str(pheno)
`data.frame':   2500 obs. of  11 variables:
 $ pca: int  1 1 1 1 1 1 1 1 1 1 ...
 $ her: int  0 1 0 0 0 0 0 0 0 0 ...
 $ age: num  67.1 70.4 64.9 60.8 64.3 ...
 $ t  : int  1 3 1 2 1 2 1 3 9 1 ...
 $ n  : int  9 9 9 9 9 9 9 0 9 9 ...
 $ m  : int  0 1 0 0 0 1 9 9 1 9 ...
 $ diffgrd: int  9 9 9 2 9 9 9 9 3 9 ...
 $ gs : int  6 7 6 7 7 7 5 6 NA 6 ...
 $ psa: num  13 75 11.3 13 9.1 51 4.3 10.7 NA 93 ...
 $ geo: int  2 2 1 1 1 1 1 1 1 1 ...
 $ ageg   : int  6 7 5 5 5 5 3 5 4 5 ...
 str(geno)
`data.frame':   2500 obs. of  125000 variables:
^C

 
 dim(geno)
[1]   2500 125000
 dim(haplo)
[1]   2500 25
 version()
Error: attempt to apply non-function
 version  
 _   
platform sparc-sun-solaris2.9
arch sparc   
os   solaris2.9  
system   sparc, solaris2.9   
status   Patched 
major2   
minor1.1 
year 2005
month07  
day  09  
language R   


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


Re: [R] Memory limits using read.table on Windows XP Pro

2005-07-19 Thread Latchezar Dimitrov
Really sorry for the wrong addressing. It was intended to the list only.
I apologize.
Latchezar 

 -Original Message-
 From: Latchezar Dimitrov 
 Sent: Tuesday, July 19, 2005 9:47 PM
 To: Latchezar Dimitrov; 'Prof Brian Ripley'
 Cc: 'r-help@stat.math.ethz.ch'
 Subject: RE: [R] Memory limits using read.table on Windows XP Pro
 
 Hello everyone,
 
 Would you please somebody explain me what my sin is (please 
 see the code and timing bellow)? And how to improve myself 
 and the following piece of R code? BTW, the code works.
 
 This is R 64-bit built by myself on Sun SPARC Solaris 9 with 
 gcc-4.0.1 (64-bit) also built by yours truly. The machine is 
 Sun Fire-V880 with 32GB memory and 8 cpu's. Nobody else was using it.
 
 
 Thank you very much
 
 Latchezar Dimitrov
 
 PS. Prof. Brian D. Ripley, thank you very much for not 
 wasting your invaluable time responding to the above. I know 
 your answer anyway. Sorry for making you read it though ...
 
   system.time(
 + haplo[i,2*j-1]-substr(as.character(geno[i,j]),1,1)
 +  ,TRUE)
 [1] 66.27 13.67 80.02  0.00  0.00
  ls()
 [1] P geno  haplo i j lrmodel pheno
  i
 [1] 1
  j
 [1] 1
  dim(P)
 [1]  1 125000
  str(P)
  num [1, 1:125000] 0.6188 0.0533 0.0893 0.8994 0.0316 ...
  str(pheno)
 `data.frame':   2500 obs. of  11 variables:
  $ pca: int  1 1 1 1 1 1 1 1 1 1 ...
  $ her: int  0 1 0 0 0 0 0 0 0 0 ...
  $ age: num  67.1 70.4 64.9 60.8 64.3 ...
  $ t  : int  1 3 1 2 1 2 1 3 9 1 ...
  $ n  : int  9 9 9 9 9 9 9 0 9 9 ...
  $ m  : int  0 1 0 0 0 1 9 9 1 9 ...
  $ diffgrd: int  9 9 9 2 9 9 9 9 3 9 ...
  $ gs : int  6 7 6 7 7 7 5 6 NA 6 ...
  $ psa: num  13 75 11.3 13 9.1 51 4.3 10.7 NA 93 ...
  $ geo: int  2 2 1 1 1 1 1 1 1 1 ...
  $ ageg   : int  6 7 5 5 5 5 3 5 4 5 ...
  str(geno)
 `data.frame':   2500 obs. of  125000 variables:
 ^C
 
  
  dim(geno)
 [1]   2500 125000
  dim(haplo)
 [1]   2500 25
  version()
 Error: attempt to apply non-function
  version
  _   
 platform sparc-sun-solaris2.9
 arch sparc   
 os   solaris2.9  
 system   sparc, solaris2.9   
 status   Patched 
 major2   
 minor1.1 
 year 2005
 month07  
 day  09  
 language R   
 

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


Re: [R] Problems with date-format (R 2.1.1 + chron)

2005-07-19 Thread Gabor Grothendieck
I think its likely that you are using different versions of chron.
I noticed that version 2.2-33 of chron had the statement:

   tms - dts - trunc(dts)

but version 2.2-35 seems to have replaced it with:

  tms - dts - floor(dts)

and that seems to be causing the problem.

As a workaround:

floor.Date - floor.trunc
days(as.Date(21-07-2005, %d-%m-%y))


On 7/19/05, Carsten Steinhoff [EMAIL PROTECTED] wrote:
 Hello,
 
 today I've updated on the newest R-Version. But sadly a function I needed
 didnt want to work:
 The input is e.g.
 
 days(as.Date(21-07-2005,%d-%m-%y))
 
 the error is: Fehler in Math.Date(dts): floor nicht definiert für Date
 Objekte
 (Error in Math.Date(dts): floor not defined for date objects)
 
 Same for year. Only months gives me the correct output.
 In Version 2.01 it worked very well, with the same chron library.
 Whats wrong ?
 
 Carsten
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


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


Re: [R] Taking the derivative of a quadratic B-spline

2005-07-19 Thread James McDermott
Would the unique quadratic defined by the three points be the same
curve as the curve predicted by a quadratic B-spline (fit to all of
the data) through those same three points?

Jim

On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
 On 7/19/2005 3:34 PM, James McDermott wrote:
  I wish it were that simple (perhaps it is and I am just not seeing
  it).  The output from cobs( ) includes the B-spline coefficients and
  the knots.  These coefficients are not the same as the a, b, and c
  coefficients in a quadratic polynomial.  Rather, they are the
  coefficients of the quadratic B-spline representation of the fitted
  curve.  I need to evaluate a linear combination of basis functions and
  it is not clear to me how to accomplish this easily.  I was hoping to
  find an alternative way of getting the derivatives.
 
 I don't know COBS, but doesn't predict just evaluate the B-spline?  The
 point of what I posted is that the particular basis doesn't matter if
 you can evaluate the quadratic at 3 points.
 
 Duncan Murdoch
 
 
  Jim McDermott
 
  On 7/19/05, Duncan Murdoch [EMAIL PROTECTED] wrote:
  On 7/19/2005 2:53 PM, James McDermott wrote:
   Hello,
  
   I have been trying to take the derivative of a quadratic B-spline
   obtained by using the COBS library.  What I would like to do is
   similar to what one can do by using
  
   fit-smooth.spline(cdf)
   xx-seq(-10,10,.1)
   predict(fit, xx, deriv = 1)
  
   The goal is to fit the spline to data that is approximating a
   cumulative distribution function (e.g. in my example, cdf is a
   2-column matrix with x values in column 1 and the estimate of the cdf
   evaluated at x in column 2) and then take the first derivative over a
   range of values to get density estimates.
  
   The reason I don't want to use smooth.spline is that there is no way
   to impose constraints (e.g. =0, =1, and monotonicity) as there is
   with COBS.  However, since COBS doesn't have the 'deriv =' option, the
   only way I can think of doing it with COBS is to evaluate the
   derivatives numerically.
 
  Numerical estimates of the derivatives of a quadratic should be easy to
  obtain accurately.  For example, if the quadratic ax^2 + bx + c is
  defined on [-1, 1], then the derivative 2ax + b, has 2a = f(1) - f(0) +
  f(-1), and b = (f(1) - f(-1))/2.
 
  You should be able to generalize this to the case where the spline is
  quadratic between knots k1 and k2 pretty easily.
 
  Duncan Murdoch
 
 


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


[R] nls

2005-07-19 Thread ekhous
Dear R-helpers,

I am trying to estimate a model that I am proposing, which consists of putting
an extra hidden layer in the Markov switching models. In the simplest case the
S(t) - Markov states - and w(t) - the extra hidden variables - are independent,
and w(t) is constant. Formally the model looks like this:
y(t)=c(1,y[t-1])%*%beta0*w+c(1,y[t-1])%*%beta1*(1-w). So I ran some simulations
to obtain the y's, and I am putting it into the nls:

res-nls(y~(a+b*x)*w+(c+d*x)*(1-w),start=list(a=1,b=0.3,c=-1,d=-0.2,w=0.5))

and the starting parameter values are similar to the ones I used for
simulations, however I am getting

Error in nlsModel(formula, mf, start) : singular gradient matrix at initial
parameter estimates

What am I doing wrong? I tried many different parameter values to no avail.
Thank you so much in advance,
Sincerely,

Eugene
McGill University

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


Re: [R] nls

2005-07-19 Thread Gabor Grothendieck
On 7/19/05, [EMAIL PROTECTED] [EMAIL PROTECTED] wrote:
 Dear R-helpers,
 
 I am trying to estimate a model that I am proposing, which consists of putting
 an extra hidden layer in the Markov switching models. In the simplest case the
 S(t) - Markov states - and w(t) - the extra hidden variables - are 
 independent,
 and w(t) is constant. Formally the model looks like this:
 y(t)=c(1,y[t-1])%*%beta0*w+c(1,y[t-1])%*%beta1*(1-w). So I ran some 
 simulations
 to obtain the y's, and I am putting it into the nls:
 
 res-nls(y~(a+b*x)*w+(c+d*x)*(1-w),start=list(a=1,b=0.3,c=-1,d=-0.2,w=0.5))
 
 and the starting parameter values are similar to the ones I used for
 simulations, however I am getting
 
 Error in nlsModel(formula, mf, start) : singular gradient matrix at initial
 parameter estimates
 

Your model is not identifiable.  You are using 5 parameters to describe 
a two dimensional model -- in fact, y is linear in x so anything beyond 
the intercept and slope are redundant, viz. a singular gradient.

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


Re: [R] Michaelis-menten equation

2005-07-19 Thread Chun-Ying Lee
Hi,

   We used known Vm and Km to simulate the data set (time, Cp) without
adding random error in there.  Yes, the line looks like very close
to a straight line.  But why can't we obtain the correct values with
fitting process?  We used optim first and then followed by using nls
to fit the model.  Thanks.

regards,
---Chun-ying Lee
 
 Hmm, sorry, no. I'm talking through a hole in my head there.
 
 Vm*y/(Km+y) makes OK sense. Vm is what you get for large y - passing
 from 1st order to 0th order kinetics. However, looking at the data
 
  plot(PKindex)
  abline(lm(conc~time,data=PKindex))
 
 shows that they are pretty much on a straight line, i.e. you are 
 in the domain of 0-order kinetics. So why are you expecting the rate
 of decrease to have changed by roughly 3/4 (from 2/3*Vm/Vd at y=2*Km
 to 1/2*Vm/Vd at y=Km when you reach 4.67)??


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

Re: [R] Problems with date-format (R 2.1.1 + chron)

2005-07-19 Thread Gabor Grothendieck
On 7/19/05, Gabor Grothendieck [EMAIL PROTECTED] wrote:
 I think its likely that you are using different versions of chron.
 I noticed that version 2.2-33 of chron had the statement:
 
   tms - dts - trunc(dts)
 
 but version 2.2-35 seems to have replaced it with:
 
  tms - dts - floor(dts)
 
 and that seems to be causing the problem.
 
 As a workaround:
 
floor.Date - floor.trunc

That should have been:

  floor.Date - trunc.Date

days(as.Date(21-07-2005, %d-%m-%y))
 
 
 On 7/19/05, Carsten Steinhoff [EMAIL PROTECTED] wrote:
  Hello,
 
  today I've updated on the newest R-Version. But sadly a function I needed
  didnt want to work:
  The input is e.g.
 
  days(as.Date(21-07-2005,%d-%m-%y))
 
  the error is: Fehler in Math.Date(dts): floor nicht definiert für Date
  Objekte
  (Error in Math.Date(dts): floor not defined for date objects)
 
  Same for year. Only months gives me the correct output.
  In Version 2.01 it worked very well, with the same chron library.
  Whats wrong ?
 
  Carsten
 
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide! 
  http://www.R-project.org/posting-guide.html
 


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


Re: [R] Regression lines for differently-sized groups on the same plot

2005-07-19 Thread Sundar Dorai-Raj


Laura M Marx wrote:
 Hi there,
   I've looked through the very helpful advice about adding fitted lines to 
 plots in the r-help archive, and can't find a post where someone has offered 
 a solution for my specific problem.  I need to plot logistic regression fits 
 from three differently-sized data subsets on a plot of the entire dataset.  
 A description and code are below:
   I have an unbalanced dataset consisting of three different species (hem, 
 yb, and sm), with unequal numbers of wood pieces in each species group.  I 
 am trying to generate a plot that will show the size of the wood piece on 
 the X axis, the probability of it having tree seedlings growing on it on the 
 Y (a binomial yes or no variable), and three fitted curves showing how the 
 probability of having tree seedlings changes with increasing wood piece size 
 for each species.
   I have no problem generating fits using GLM, and no problem creating the 
 plot.  However, if I try to add a fitted curve based only on the hem data 
 subset to a plot that shows the entire dataset, I get an error message that 
 the lengths of those data sets differ. Error in xy.coords(x,y) : x and y 
 lengths differ.  I could see R's point -- you can't plot a regression line 
 of babies born as a function of stork abundance on a graph of cherries 
 produced (Y) versus rainfall (X), which for all the program knows, I'm 
 trying to do.  As a temporary fix, I added NAs to the end of the hem, yb, 
 and sm subsets to make them the same length as the entire dataset.  I can 
 now add my fitted curves to the plot, but the lines are not connected.  That 
 is, if the hem group only contains wood pieces that are 1, 4, and 10 meters 
 long, the plot has an X axis that ranges from 1 to 10, but line segments for 
 the hem group regression line only appear above 1, 4, and 10.  How can I fix 
 this?  An ideal solution would not require me to make the hem subset of my 
 data the same length as the full dataset, either (although the summaries of 
 regressions with the NAs (or zeroes) added and taken away are identical).  
 I'd also settle for a work-around that would have R connect the pieces of 
 the curve so that I get a solid line rather than small dots and dashes where 
 actual data exist.  Thanks so much for your help!
   Laura Marx
   Michigan State University, Dept. of Forestry 
 
 #Note: hemdata has all the rows that are not hemlock species replaced with 
 #NAs.
 hemhem=glm(hempresence~logarea, family=binomial(logit), data=hemdata)
 hemyb=glm(hempresence~logarea, family=binomial(logit), data=birchdata)
 hemsm=glm(hempresence~logarea, family=binomial(logit), data=mapledata) 
 
 attach(logreg) #logreg is the full dataset
 plot(logarea, hempresence, xlab = Surface area of log (m2), 
 ylab=Probability of hemlock seedling presence, type=n, font.lab=2, 
 cex.lab=1.5, axes=TRUE)
 lines(logarea,fitted(hemhem), lty=1, lwd=2)
 lines(logarea,fitted(hemyb), lty=dashed, lwd=2)
 lines(logarea,fitted(hemsm), lty=dotted, lwd=2)
 

Hi, Laura,

Would ?predict.glm be better?

plot(logarea, hempresence,
  xlab = Surface area of log (m2),
  ylab=Probability of hemlock seedling presence,
  type=n, font.lab=2, cex.lab=1.5, axes=TRUE)
lines(logarea, predict(hemhem, logreg, response), lty=1, lwd=2)
lines(logarea, predict(hemyb, logreg, response), lty=dashed, lwd=2)
lines(logarea, predict(hemsm, logreg, response), lty=dotted, lwd=2)

Without seeing more description of your data, this is still a guess.

--sundar

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


[R] Turning off return warning messages.

2005-07-19 Thread Steve Su
Dear All,

Is there a way I can turn off the following warning message for using 
multi-argument returns?

multi-argument returns are deprecated in: return(p1, p2, p3, p4)

Steve.

**

 Steve Su ([EMAIL PROTECTED])
 Postdoctoral fellow

 Faculty of Business 
 Queensland University of Technology  
 
 Postal Address: Steve Su, School of Accountancy, QUT, PO Box 2434, Brisbane, 
 Queensland, Australia, 4001.   

 CRICOS No. 00213J 
 
 Phone:  +61 7 3864 4357
 Fax:+61 7 3864 1812  
 Mobile: 0421  840  586  
 .   
   _--_|\  
  /  QUT   
  \_.--._/  
 
v   
  
**
   
[[alternative HTML version deleted]]

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


[R] Question about Installing SJava package

2005-07-19 Thread Taemyong Choi
Dear all,

I have an error message installing SJava package.
So I searched web site(google) and R-mailing list to find a similar error 
message.
But I couldn't find it.

I installed R-2.1.1 like this on Fedora Core4

1) /configure --enable-R-shlib --with-libpng --with-jpeglib
2) make - make check - make install

and then issuing on shell prompt (red lines are error messages)

R CMD INSTALL -c /usr/local/src/R/SJava_0.68-0.tar.gz

* Installing *source* package 'SJava' ...
checking for java... /usr/java/jdk1.5.0_04//bin/java
Java VM /usr/java/jdk1.5.0_04//bin/java
checking for javah... /usr/java/jdk1.5.0_04//bin/javah
Looking in /usr/java/jdk1.5.0_04/include
Looking in /usr/java/jdk1.5.0_04/include/linux
checking for g++... g++
checking for C++ compiler default output... a.out
checking whether the C++ compiler works... yes
checking whether we are cross compiling... no
checking for suffix of executables...
checking for suffix of object files... o
checking whether we are using the GNU C++ compiler... yes
checking whether g++ accepts -g... yes
checking for gcc... gcc
checking whether we are using the GNU C compiler... yes
checking whether gcc accepts -g... yes
checking for gcc option to accept ANSI C... none needed
checking for Rf_initEmbeddedR in -lR... no
No R shared library found
configure: creating ./config.status
config.status: creating Makevars
config.status: creating src/Makevars
config.status: creating src/RSJava/Makefile
config.status: creating Makefile_rules
config.status: creating inst/scripts/RJava.bsh
config.status: creating inst/scripts/RJava.csh
config.status: creating R/zzz.R
config.status: creating cleanup
config.status: creating inst/scripts/RJava
Copying the cleanup script to the scripts/ directory
Building libRSNativeJava.so in /tmp/R.INSTALL.tf2988/SJava/src/RSJava
if test ! -d /usr/local/lib/R/library/SJava/libs ; then \
mkdir /usr/local/lib/R/library/SJava/libs ; \
fi
gcc -g -O2 -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-c CtoJava.c
CtoJava.cweb:215: error: static declaration of 'std_env' follows
non-static declaration
CtoJava.cweb:195: error: previous declaration of 'std_env' was here
make: *** [CtoJava.o] Error 1
Generating JNI header files from Java classes.
RForeignReference, RManualFunctionActionListener, ROmegahatInterpreter 
REvaluator
*
Warning:
At present, to use the library you must set the
LD_LIBRARY_PATH environment variable
to
/usr/local/lib/R/library/SJava/libs:/usr/java/jdk1.5.0_04/jre/lib/i386/client:/usr/java/jdk1.5.0_04/jre/lib/i386:/usr/java/jdk1.5.0_04/jre/../lib/i386:
or use one of the RJava.bsh or RJava.csh scripts
*
** libs
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include -IRSJava
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-I/usr/local/include -fPIC -g -O2 -c ConverterExamples.c -o
ConverterExamples.o
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include -IRSJava
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-I/usr/local/include -fPIC -g -O2 -c Converters.c -o Converters.o
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include -IRSJava
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-I/usr/local/include -fPIC -g -O2 -c REmbed.c -o REmbed.o
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include -IRSJava
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-I/usr/local/include -fPIC -g -O2 -c REmbedWin.c -o REmbedWin.o
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include -IRSJava
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-I/usr/local/include -fPIC -g -O2 -c REval.c -o REval.o
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.
-I/tmp/R.INSTALL.tf2988/SJava/inst/include -IRSJava
-I/usr/java/jdk1.5.0_04//include -I/usr/java/jdk1.5.0_04//include/linux
-I/usr/local/include -fPIC -g -O2 -c RFunctionListener.c -o
RFunctionListener.o
gcc -I/usr/local/lib/R/include -D_R_ -I/usr/local/lib/R/include
-I/usr/local/lib/R/include/R_ext
-I/tmp/R.INSTALL.tf2988/SJava/src/RSJava -I.

[R] class in apply

2005-07-19 Thread Omar Lakkis
Numeric data that is part of a mixed type data frame is converted into
character. How can I tell apply to maintain the original class of a
column and not convert it into character. I would like to do this of
the vector and not inside the apply function individually over each
element. Consider the following two scenarios, in the second column
'a' maintained its class while it lost its numeric type in the first
case.

 df = data.frame(a=c(1,2), b=c('A','B'))
 df
  a b
1 1 A
2 2 B
 a=apply(df, 1, function(r) print(class(r['a'])))
[1] character
[1] character
 a=apply(df, 1, function(r) print(class(r['b'])))
[1] character
[1] character


 df = data.frame(a=c(1,2))
 df
  a
1 1
2 2
 a=apply(df, 1, function(r) print(class(r['a'])))
[1] numeric
[1] numeric

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


Re: [R] Turning off return warning messages.

2005-07-19 Thread Paul Roebuck
On Wed, 20 Jul 2005, Steve Su wrote:

 Is there a way I can turn off the following warning message
 for using multi-argument returns?

 multi-argument returns are deprecated in: return(p1, p2, p3, p4)

doubleEm - function(p1, p2, p3, p4) {
return(list(p1 = p1*p1,
p2 = p2*p2,
p3 = p3*p3,
p4 = p4*p4))
}

res - doubleEm(1, 2, 3, 4)
cat(p3 =, res$p3, \n)

--
SIGSIG -- signature too long (core dumped)

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