Re: [R] How to convert c:\a\b to c:/a/b

2005-06-29 Thread Prof Brian Ripley
On Wed, 29 Jun 2005, David Duffy wrote:

 I couldn't resist adding a more literal answer

This can only work for escapes which are preserved.  The parser maps
\n to a character (LF) and the deparser maps it back to \n.
This happens to be true of \a \b \f \n \r \t \v \\ but no others.

For example, \s is mapped to s, and there is no difference between \s and
s in the parsed input.


 unback - function(x) {
  chars - unlist(strsplit(deparse(x),))
  chars - chars[-c(1,length(chars))]
  paste(gsub(,/,chars),collapse=)
 }

 unback(\n)

 unback(\s)
[1] s

Spencer Graves keeps on insisting there is a better way, but all the
solutions are to avoid sending the string to the parser, and hence
avoiding having the string directly in an R script.  This is common in 
shell scripts, which use 'here' documents to avoid 'quoting hell'.

We can do that in R too. Here are two variants I have not seen in the 
thread

test1.R:
scan(, , allowEscapes=FALSE, n=1, quiet=TRUE)
D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
catIx, \n, sep=)

R --slave --vanilla  test1.R
D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt

(This one does not allow quoted strings.)

test2.R:
x - readLines(stdin(), n=1)
D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
x - gsub('^(.*)$', \\1, x)
cat(x, \n)

R --slave --vanilla  test2.R
D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt

(This one allows surrounding quotes or not.)

-- 
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] where can i download the metrics package?

2005-06-29 Thread Laura Holt
There are several packages within rmetrics such as fSeries, fBasics, 
fExtremes, and so on.

You can download those in the usual way.


From: Spencer Graves [EMAIL PROTECTED]
To: ronggui [EMAIL PROTECTED]
CC: R [EMAIL PROTECTED]
Subject: Re: [R] where can i download the metrics package?
Date: Tue, 28 Jun 2005 08:27:47 -0700

 How about Rmetrics, listed under misc:  related projects on
www.r-project.org.

 spencer graves

ronggui wrote:

  i learn it from metrics: Towards a package for doing econometrics in 
Rbut i can not find it in cran.
 
 

--
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-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] quick way to construct formula

2005-06-29 Thread Eric Lecoutre


Here is a way: it uses 'paste' but I dont think it is a problem anyway
to use it.
Nevertheless, it is surely a bad idea to fit any model with more than
25 terms...


   main_effects = paste(nam,collapse=+)
   inter - outer(nam,nam,paste,sep=:)
   inter - paste(inter[upper.tri(inter)],collapse=+)
   log_effects - paste(log(,nam,),sep=,collapse=+)

as.formula(paste(~,main_effects,+,inter,+,log_effects,sep=))
~x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x1:x2 + x1:x3 + 
x2:x3 + x1:x4 + x2:x4 + x3:x4 + x1:x5 + x2:x5 + x3:x5 + x4:x5 + 
x1:x6 + x2:x6 + x3:x6 + x4:x6 + x5:x6 + x1:x7 + x2:x7 + x3:x7 + 
x4:x7 + x5:x7 + x6:x7 + x1:x8 + x2:x8 + x3:x8 + x4:x8 + x5:x8 + 
x6:x8 + x7:x8 + x1:x9 + x2:x9 + x3:x9 + x4:x9 + x5:x9 + x6:x9 + 
x7:x9 + x8:x9 + x1:x10 + x2:x10 + x3:x10 + x4:x10 + x5:x10 + 
x6:x10 + x7:x10 + x8:x10 + x9:x10 + log(x1) + log(x2) + log(x3) + 
log(x4) + log(x5) + log(x6) + log(x7) + log(x8) + log(x9) + 
log(x10)




HTH,

Eric



Eric Lecoutre
UCL /  Institut de Statistique
Voie du Roman Pays, 20
1348 Louvain-la-Neuve
Belgium

tel: (+32)(0)10473050
[EMAIL PROTECTED]
http://www.stat.ucl.ac.be/ISpersonnel/lecoutre

If the statistics are boring, then you've got the wrong numbers. -Edward
Tufte   


 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of Luke
 Sent: mercredi 29 juin 2005 4:49
 To: R-help@stat.math.ethz.ch
 Subject: [R] quick way to construct formula
 
 
 Dear R users,
 
 I have a data with 1000 variables named x1, x2, ..., 
 x1000, and I want to construct a formula like this format:
 
 ~x1+x2+...+x1000+x1:x2+x1:x3+x999:x1000+log(x1)+...+log(x1000)
 
 That is: the base variables followed by all interaction terms 
 and all base feature log-transformations. I know I can use 
 several paste functions to construct it. But is there any 
 other handy way to do it?
 
 -Luke
 
 __
 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] enhanced MDS

2005-06-29 Thread falissard
Try also ?sphpca (library(psy) -- new version 0.7)
Best
Bruno



Bruno Falissard
INSERM U669, PSIGIAM
Paris Sud Innovation Group in Adolescent Mental Health
Maison de Solenn
97 Boulevard de Port Royal
75679 Paris cedex 14, France
tel : (+33) 6 81 82 70 76
fax : (+33) 1 45 59 34 18
web site : http://perso.wanadoo.fr/bruno.falissard/

 

-Message d'origine-
De : [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] De la part de Karen Kotschy
Envoyé : mardi 28 juin 2005 13:39
À : r-help@stat.math.ethz.ch
Objet : [R] enhanced MDS

Hi again

Sorry, in looking again at sammon and isoMDS I see that they seem to do 
exactly what I want, except that they are non-metric, which means, as I 
understand it, that they relate the rank orders of the variables rather than

the actual distances. 

Could I use these non-metric MDS packages even if my distances are metric?

Thanks
Karen
-- 
Karen Kotschy
Centre for Water in the Environment
University of the Witwatersrand
Johannesburg

P/Bag X3, Wits, 2050
Tel: +2711 717-6425

__
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] quick way to construct formula

2005-06-29 Thread Prof Brian Ripley
On Tue, 28 Jun 2005, Luke wrote:

 Dear R users,

 I have a data with 1000 variables named x1, x2, ..., x1000, and
 I want to construct a formula like this format:

 ~x1+x2+...+x1000+x1:x2+x1:x3+x999:x1000+log(x1)+...+log(x1000)

 That is: the base variables followed by all interaction terms and all
 base feature log-transformations. I know I can use several paste
 functions to construct it. But is there any other handy way to do it?

Do you really want a formula with over half a million terms in?  I think 
it is likely that R will hit recursion limits in handing such a formula, 
quite possibly C-level stack overflow.

For a more modest example:

dd - data.frame(x1=1, x2=2, x3=3)
terms(~ . + .^2, data=dd, simplify=TRUE)

will get you all terms and their interactions.  Using paste() to get the 
log terms is the easiest way I know.


-- 
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] enhanced multidimensional scaling?

2005-06-29 Thread Martin Maechler
 Karen == Karen Kotschy [EMAIL PROTECTED]
 on Tue, 28 Jun 2005 13:13:31 +0200 writes:

Karen Dear R list Would anyone be able to tell me whether
Karen it is possible to do enhanced multidimensional
Karen scaling (enhanced MDS) in R? In other words,
Karen something that goes beyond cmdscale by iteratively
Karen improving the fit between observed dissimilarities
Karen and inter-object distances, using the KYST algorithm
Karen (Kruskal, 1964).

Since you know about cmdscale(), do read it's help page more
carefully, in particular the section  See Also:
{which, in an ESS help buffer, is reachable by simple s s
The help page of the first function mentioned there is entitled
`` Kruskal's Non-metric Multidimensional Scaling ''

Karen I have found several implementations of non-metric
Karen MDS in various packages but nothing like what I have
Karen described above.

well, you've mentioned Kruskal(1964), but not described much
more, hence I hope the above hint does help.

Regards,
Martin Maechler, ETH Zurich

__
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] comparison of packages for Unit Root test

2005-06-29 Thread Amir Safari
 
 
Dear R Users,
 
Could somebody please compare the packages of unit root test ( Uroot, Ucra, 
tseries and fseries ) regarding the type of test ( without constant and trend, 
with constant , and with constant and trend ) ?
 
Regards,
Amir Safari
 


-

 Rekindle the Rivalries. Sign up for Fantasy Football
[[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] Running SVM {e1071}

2005-06-29 Thread Amir Safari
 
Dear David, Dear Friends,
 
After any running svm I receive different results of Error estimation of 'svm' 
using 10-fold cross validation. What is the reason ? It is caused by the 
algorithm, libsvm , e1071 or something els? Which value can be optimal one ? 
How much run can reach to the optimality.And finally, what is difference 
between Error estimation of svm using 10-fold cross validation and MSE ( Mean 
Square Error ) ?
So many thanks in advance for your help.
Best Regards,
Amir 
 
 

__



[[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] comparison of packages for Unit Root test

2005-06-29 Thread Achim Zeileis
On Wed, 29 Jun 2005 03:08:23 -0700 (PDT) Amir Safari wrote:

  
  
 Dear R Users,
  
 Could somebody please compare the packages of unit root test ( Uroot,
 Ucra, tseries and fseries ) regarding the type of test ( without
 constant and trend, with constant , and with constant and trend ) ?

IMHO, the tseries implementation is still the reference implementation
for the significance tests which also have a slim interface.
urca (sic!) offers some more tools in an integrated way. If you are not
just interested in one test or the other, urca certainly offers much
tools in a unit root analysis.
If you use the Rmetrics suite of packages, then fSeries (sic!) is
probably most useful.
uroot (sic!) offers a GUI which is a bit cumbersome to set up and the
output is a bit crude, e.g., doesn't offer p-values.

my 2 cent,
Z

 Regards,
 Amir Safari
  
 
   
 -
 
  Rekindle the Rivalries. Sign up for Fantasy Football
   [[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


Re: [R] Running SVM {e1071}

2005-06-29 Thread Uwe Ligges
Amir Safari wrote:
  
 Dear David, Dear Friends,
  
 After any running svm I receive different results of Error estimation of 
 'svm' using 10-fold cross validation. What is the reason ? It is caused by 
 the algorithm, libsvm , e1071 or something els? Which value can be optimal 
 one ? How much run can reach to the optimality.And finally, what is 
 difference between Error estimation of svm using 10-fold cross validation and 
 MSE ( Mean Square Error ) ?
 So many thanks in advance for your help.

10-fold cross validation chooses the traing/test sets randomly, hence 
you won't get the same results each time you run it ...

Uwe Ligges


 Best Regards,
 Amir 
  
  
 
 __
 
 
 
   [[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


Re: [R] svm and scaling input

2005-06-29 Thread David Meyer

[EMAIL PROTECTED] wrote:
 Dear All,
 
 I've a question about scaling the input variables for an analysis with
 svm (package e1071). Most of my variables are factors with 4 to 6
 levels but there are also some numeric variables.
 
 I'm not familiar with the math behind svms, so my assumtions maybe
 completely wrong ... or obvious. Will the svm automatically expand the
 factors into a binary matrix? 

yes.

 If I add numeric variables outside the range of 0 to 1 do I have to
 scale them to have 0 to 1 range?

svm() will scale your data by default.

Cheers, 
David

-- 
Dr. David Meyer
Department of Information Systems and Operations

Vienna University of Economics and Business Administration
Augasse 2-6, A-1090 Wien, Austria, Europe
Fax: +43-1-313 36x746 
Tel: +43-1-313 36x4393
HP:  http://wi.wu-wien.ac.at/~meyer/

__
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] Running SVM {e1071}

2005-06-29 Thread David Meyer

  
 Dear David, Dear Friends,
  
 After any running svm I receive different results of Error estimation
 of 'svm' using 10-fold cross validation. 

using tune.svm(), or the `cross' parameter of svm()?

 What is the reason ? It is caused by the algorithm, libsvm , e1071 or
 something els? 

The splits are chosen randomly.

 Which value can be optimal one? 

The Bayes Error.

 How much run can reach to the optimality.

What do you mean by `How much run'?

 And finally, what is difference between Error estimation of svm using
 10-fold cross validation and MSE ( Mean Square Error ) ?

the former is an error estimation _procedure_, the latter is an error
_measure.

Cheers,
David

__
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] Memory Management under Linux: Problems to allocate large amounts of data

2005-06-29 Thread Dubravko Dolic
Dear Group

I'm still trying to bring many data into R (see older postings). After solving 
some troubles with the database I do most of the work in MySQL. But still I 
could be nice to work on some data using R. Therefore I can use a dedicated 
Server with Gentoo Linux as OS hosting only R. This Server is a nice machine 
with two CPU and 4GB RAM which should do the job:

Dual Intel XEON 3.06 GHz
4 x 1 GB RAM PC2100 CL2
HP Proliant DL380-G3

I read the R-Online help on memory issues and the article on garbage collection 
from the R-News 01-2001 (Luke Tierney). Also the FAQ and some newsgroup 
postings were very helpful on understanding memory issues using R.

Now I try to read data from a database. The data I wanted to read consists of 
158902553 rows and one field (column) and is of type bigint(20) in the 
database. I received the message that R could not allocate the 2048000 Kb 
(almost 2GB) sized vector. As I have 4BG of RAM I could not imagine why this 
happened. In my understanding R under Linux (32bit) should be able to use the 
full RAM. As there is not much space used by OS and R as such (free shows the 
use of app. 670 MB after dbSendQuery and fetch) there are 3GB to be occupied by 
R. Is that correct?

After that I started R by setting n/vsize explicitly

R --min-vsize=10M --max-vsize=3G --min-nsize=500k --max-nsize=100M

 mem.limits()
nsize vsize
104857600NA

and received the same message.


A garbage collection delivered the following information:

 gc()
 used (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
Ncells 217234  5.9 50   13.4   280050   13.4
Vcells  87472  0.7  157650064 1202.8   3072 196695437 1500.7


Now I'm at a loss. Maybe anyone could give me a hint where I should read 
further or which Information can take me any further 





Dubravko Dolic
Statistical Analyst
Tel:      +49 (0)89-55 27 44 - 4630
Fax:     +49 (0)89-55 27 44 - 2463
Email: [EMAIL PROTECTED]
Komdat GmbH
Nymphenburger Straße 86
80636 München
-
ONLINE MARKETING THAT WORKS
-
This electronic message contains information from Komdat Gmb...{{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] coxph, survfit and Brier score

2005-06-29 Thread Stephen Henderson
Hello and apologies for a very long question. I thought it better to be
verbose and clear than short and imprecise. I am trying to compute the brier
score comparing true surv object of test data to predictions from train data
(using sbrier in ipred package). I am having trouble getting the right
format for my prediction object I think. I have split the DLBCL dataset and
fit a coxph on a training set then I am trying to get out a prediction using
the survfit function with newdata=test. My code is below. 

library(survival)
library(ipred)


data(DLBCL)
DLBCL[complete.cases(DLBCL),]-DLBCL
train-DLBCL[1:19,]
test-DLBCL[20:38,]
train.surv-Surv(train$time, train$cens)
test.surv-Surv(test$time, test$cens)


train.mod-coxph(train.surv~IPI, data=train)
 pred- survfit(train.mod, newdata=test)

class(pred)
[1] survfit.cox survfit

pred

Call: survfit.coxph(object = train.mod, newdata = test)

   n events median 0.95LCL 0.95UCL
 [1,] 19 10Inf27.1 Inf
 [2,] 19 10Inf Inf Inf
 [3,] 19 10Inf Inf Inf
 [4,] 19 10   23.7 4.1 Inf
 [5,] 19 10   23.7 4.1 Inf
 [6,] 19 10Inf27.1 Inf
 [7,] 19 10Inf Inf Inf
 [8,] 19 10Inf Inf Inf
 [9,] 19 10Inf27.1 Inf
[10,] 19 104.1 2.9 Inf
[11,] 19 10Inf27.1 Inf
[12,] 19 104.1 2.9 Inf
[13,] 19 10Inf Inf Inf
[14,] 19 10   23.7 4.1 Inf
[15,] 19 10Inf27.1 Inf
[16,] 19 10Inf27.1 Inf
[17,] 19 10Inf Inf Inf
[18,] 19 10Inf Inf Inf
[19,] 19 10Inf27.1 Inf


sbrier(test.surv, pred)
Error in switch(ptype, survfit = { : switch: EXPR must return a length 1
vector

The sbrier function clearly does not recognize the format of this survfit
object pred. Indeed it seems to have the same variables but does not look
like the object you receive from a normal survfit call e.g. 


KM - survfit(train.surv)
KM
Call: survfit(formula = train.surv)

  n  events  median 0.95LCL 0.95UCL 
   19.010.071.315.5 Inf 

sbrier(train.surv, KM)
integrated Brier score 
 0.2220228 
attr(,time)
[1]   2.4 102.4

My question is can I either force survfit(etc..., newdata=x) to return a
useful survfit object for use with sbrier, or alternatively coerce the pred
object above into a survfit object? Or am I missing something ..is there a
reason survfit returns the predicted object in a different format?

PS I realize the fit is not good but just want to get the code to work.

Thank You

Stephen Henderson


**
This email and any files transmitted with it are confidentia...{{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


Re: [R] Producing character given i.e. | with plotmath

2005-06-29 Thread Gorjanc Gregor
 Hello!
 
 Does someone know how to produce
 
   L(y|mu)
 
 with plotmath?
 
 Some code with unsuccessfull results:
 
 plot(dnorm(x = seq(from = -4, to = 4, by = 0.1)), type = l)
 ## Not what I want
 legend(legend = c(expression(L(y:mu))), x = topright)
 
 ## Strange, is this a bug?
 legend(legend = c(expression(L(y|mu))), x = top)

 No, | is a logical Operator that can be rewritten in its original 
 function form as follows:

 |(FALSE, TRUE)

 Hence the result is expected.

Yes, that makes sense. Thanks!

 ## Group produces an error
 legend(legend = c(expression(group(L(y, |, mu, x = topleft)

 You have not specified any delimiter.
Ok, I got this point wrong from ?plotmath .

 ## Paste keeps commas in expression
 legend(legend = c(expression(paste(L(y, |, mu, x = bottomleft)

correct

 ## This one is OK, but braces are not as they should be 
 legend(legend = c(expression(paste(L(y, |, mu, x = bottom)

 What's wrong with the braces?`
They are not the same as in previous cases. They are not bold.

 What you really want is:
legend(legend = c(expression(L(group(, y, |) * mu))),
  x = center)

Yes, that is what I want. Thanks!

I got additionally response from Roger D. Peng, which stated:

 How about
 legend(legend = expression(L(y |  * mu * )), x = topleft)

That's also OK, but note above comment about bold braces.

Regards, Gregor

__
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] x*x*x*... vs x^n

2005-06-29 Thread Robin Hankin
Hi

I have been wondering if there one can speed up calculating small powers
of numbers such as x^8 using multiplication.

In addition, one can be a bit clever and calculate x^8 using only 3  
multiplies.

look at this:


  f1 - function(x){x*x*x*x*x*x*x*x}
  f2 - function(x){x^8}
  f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}

[so f1() and f2() and f3() are algebraically identical]


  a - rnorm(100)
  system.time(ignore - f1(a))
[1] 0.50 0.17 2.88 0.00 0.00

  system.time(ignore - f2(a))
[1] 0.31 0.03 1.40 0.00 0.00

  system.time(ignore - f3(a))
[1] 0.10 0.07 0.18 0.00 0.00


[these figures show little variance from trial to trial]


I was expecting f2() and f3() to be about the same.
I was not expecting a factor of 3 there!

anyone got any comments?




--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

__
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] x*x*x*... vs x^n

2005-06-29 Thread Ken Knoblauch
Something like this is exploited very nicely in the mtx.exp 
for matrix powers in the Malmig package, actually.



Hi

I have been wondering if there one can speed up calculating small powers
of numbers such as x^8 using multiplication.

In addition, one can be a bit clever and calculate x^8 using only 3  
multiplies.

look at this:


  f1 - function(x){x*x*x*x*x*x*x*x}
  f2 - function(x){x^8}
  f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}

[so f1() and f2() and f3() are algebraically identical]


  a - rnorm(100)
  system.time(ignore - f1(a))
[1] 0.50 0.17 2.88 0.00 0.00

  system.time(ignore - f2(a))
[1] 0.31 0.03 1.40 0.00 0.00

  system.time(ignore - f3(a))
[1] 0.10 0.07 0.18 0.00 0.00


[these figures show little variance from trial to trial]


I was expecting f2() and f3() to be about the same.
I was not expecting a factor of 3 there!

anyone got any comments?




--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743


Ken Knoblauch
Inserm U371, Cerveau et Vision
Department of Cognitive Neurosciences
18 avenue du Doyen Lepine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: 06 84 10 64 10
http://www.lyon.inserm.fr/371/

__
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] x*x*x*... vs x^n

2005-06-29 Thread Tuszynski, Jaroslaw W.
I tried your code and got different results:
 system.time(ignore - f1(a))
[1] 0.83 0.09 1.08   NA   NA
 system.time(ignore - f2(a))
[1] 0.38 0.01 0.41   NA   NA
 system.time(ignore - f3(a))
[1] 0.32 0.04 0.43   NA   NA

So I tried it again but with a loop and got:

  for(i in 1:10) cat(system.time(ignore - f2(a)), \n)
0.36 0.04 0.44 NA NA 
0.32 0.01 0.34 NA NA 
0.28 0.03 0.32 NA NA 
0.29 0.03 0.35 NA NA 
0.3 0.02 0.32 NA NA 
0.28 0.03 0.32 NA NA 
0.3 0.02 0.32 NA NA 
0.29 0.02 0.34 NA NA 
0.23 0.03 0.32 NA NA 
0.42 0 0.45 NA NA  

 for(i in 1:10) cat(system.time(ignore - f3(a)), \n)
0.19 0.04 0.25 NA NA 
0.17 0.04 0.25 NA NA 
0.21 0.02 0.25 NA NA 
0.21 0.02 0.23 NA NA 
0.18 0.04 0.23 NA NA 
0.18 0.05 0.23 NA NA 
0.18 0.04 0.25 NA NA 
0.17 0.06 0.23 NA NA 
0.2 0.02 0.23 NA NA 
0.14 0.06 0.25 NA NA 

It seems to me that f3 is 50% slower than f2 not 300%.

My System is:
- R version: R 2.1.1
- Operating System: Win XP
- Compiler: mingw32-gcc-3.4.2

Jarek
\===

 Jarek Tuszynski, PhD.   o / \ 
 Science Applications International Corporation  \__,|  
 (703) 676-4192  \
 [EMAIL PROTECTED] `\



-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Robin Hankin
Sent: Wednesday, June 29, 2005 7:32 AM
To: r-help
Subject: [R] x*x*x*... vs x^n

Hi

I have been wondering if there one can speed up calculating small powers of
numbers such as x^8 using multiplication.

In addition, one can be a bit clever and calculate x^8 using only 3
multiplies.

look at this:


  f1 - function(x){x*x*x*x*x*x*x*x}
  f2 - function(x){x^8}
  f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}

[so f1() and f2() and f3() are algebraically identical]


  a - rnorm(100)
  system.time(ignore - f1(a))
[1] 0.50 0.17 2.88 0.00 0.00

  system.time(ignore - f2(a))
[1] 0.31 0.03 1.40 0.00 0.00

  system.time(ignore - f3(a))
[1] 0.10 0.07 0.18 0.00 0.00


[these figures show little variance from trial to trial]


I was expecting f2() and f3() to be about the same.
I was not expecting a factor of 3 there!

anyone got any comments?




--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton European Way, Southampton SO14
3ZH, UK
  tel  023-8059-7743

__
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] x*x*x*... vs x^n

2005-06-29 Thread Duncan Murdoch
On 6/29/2005 7:32 AM, Robin Hankin wrote:
 Hi
 
 I have been wondering if there one can speed up calculating small powers
 of numbers such as x^8 using multiplication.
 
 In addition, one can be a bit clever and calculate x^8 using only 3  
 multiplies.
 
 look at this:
 
 
   f1 - function(x){x*x*x*x*x*x*x*x}
   f2 - function(x){x^8}
   f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}
 
 [so f1() and f2() and f3() are algebraically identical]
 
 
   a - rnorm(100)
   system.time(ignore - f1(a))
 [1] 0.50 0.17 2.88 0.00 0.00
 
   system.time(ignore - f2(a))
 [1] 0.31 0.03 1.40 0.00 0.00
 
   system.time(ignore - f3(a))
 [1] 0.10 0.07 0.18 0.00 0.00
 
 
 [these figures show little variance from trial to trial]
 
 
 I was expecting f2() and f3() to be about the same.
 I was not expecting a factor of 3 there!
 
 anyone got any comments?

If you look in src/main/arithmetic.c, you'll see that R does a general 
real-valued power (using C's pow() function) whenever either one of the 
args is real (except for a few special cases, e.g. non-numbers, or 
powers of 2 or 0.5).  There is an internal R function equivalent to your 
f3, but it is not used in the situation of real^integer (and in any 
case, x^8 is real^real).

I suppose if you wanted to submit a patch someone would take a look, but 
the question is whether there is really any calculation whose execution 
time would be materially affected by this.  Most computations are not 
dominated by integer power calculations, so is this really worth the 
trouble?

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] range of input data

2005-06-29 Thread Carsten Steinhoff
Hi,
 
I've written this function:
 
g = function(test,p1,p2)
{ 
 
test=sort(test)
merke=0
for (z in 1:length(test))
{
  F1=((2*z-1)/length(test))
  F21=log(plnorm(test[z],p1,p2))
  F22=log(1-plnorm(test[length(test)+1-z],p1,p2))
  F2= F21+F22
  merke=merke+F1*F2
 
}
return(-length(test)*-merke)
}
 
 
When I call it for e.g. test = 1:10 it only returns ONE value. Instead what
I want to have is 10 values (like it makes the dlnorm for instance).
I could make that with a loop, but maybe there is there a more simple trick
?
 
Carsten

[[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] Memory Management under Linux: Problems to allocate large amounts of data

2005-06-29 Thread Prof Brian Ripley
Let's assume this is a 32-bit Xeon and a 32-bit OS (there are 
64-bit-capable Xeons).  Then a user process like R gets a 4GB address 
space, 1GB of which is reserved for the kernel.  So R has a 3GB address 
space, and it is trying to allocate a 2GB contigous chunk.  Because of 
memory fragmentation that is quite unlikely to succeed.

We run 64-bit OSes on all our machines with 2GB or more RAM, for this 
reason.

On Wed, 29 Jun 2005, Dubravko Dolic wrote:

 Dear Group

 I'm still trying to bring many data into R (see older postings). After 
 solving some troubles with the database I do most of the work in MySQL. 
 But still I could be nice to work on some data using R. Therefore I can 
 use a dedicated Server with Gentoo Linux as OS hosting only R. This 
 Server is a nice machine with two CPU and 4GB RAM which should do the 
 job:

 Dual Intel XEON 3.06 GHz
 4 x 1 GB RAM PC2100 CL2
 HP Proliant DL380-G3

 I read the R-Online help on memory issues and the article on garbage 
 collection from the R-News 01-2001 (Luke Tierney). Also the FAQ and some 
 newsgroup postings were very helpful on understanding memory issues 
 using R.

 Now I try to read data from a database. The data I wanted to read 
 consists of 158902553 rows and one field (column) and is of type 
 bigint(20) in the database. I received the message that R could not 
 allocate the 2048000 Kb (almost 2GB) sized vector. As I have 4BG of RAM 
 I could not imagine why this happened. In my understanding R under Linux 
 (32bit) should be able to use the full RAM. As there is not much space 
 used by OS and R as such (free shows the use of app. 670 MB after 
 dbSendQuery and fetch) there are 3GB to be occupied by R. Is that 
 correct?

Not really.  The R executable code and the Ncells are already in the 
address space, and this is a virtual memory OS, so the amount of RAM is 
not relevant (it would still be a 3GB limit with 12GB of RAM).

 After that I started R by setting n/vsize explicitly

 R --min-vsize=10M --max-vsize=3G --min-nsize=500k --max-nsize=100M

 mem.limits()
nsize vsize
 104857600NA

 and received the same message.


 A garbage collection delivered the following information:

 gc()
 used (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
 Ncells 217234  5.9 50   13.4   280050   13.4
 Vcells  87472  0.7  157650064 1202.8   3072 196695437 1500.7


 Now I'm at a loss. Maybe anyone could give me a hint where I should read 
 further or which Information can take me any further

-- 
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] x*x*x*... vs x^n

2005-06-29 Thread Robin Hankin
Hi  Duncan


On Jun 29, 2005, at 02:04 pm, Duncan Murdoch wrote:

 On 6/29/2005 7:32 AM, Robin Hankin wrote:

 Hi
 I have been wondering if there one can speed up calculating small  
 powers
 of numbers such as x^8 using multiplication.
 In addition, one can be a bit clever and calculate x^8 using only  
 3  multiplies.
 look at this:
   f1 - function(x){x*x*x*x*x*x*x*x}
   f2 - function(x){x^8}
   f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}


[snip]


 If you look in src/main/arithmetic.c, you'll see that R does a  
 general real-valued power (using C's pow() function) whenever  
 either one of the args is real (except for a few special cases,  
 e.g. non-numbers, or powers of 2 or 0.5).  There is an internal R  
 function equivalent to your f3, but it is not used in the situation  
 of real^integer (and in any case, x^8 is real^real).

 I suppose if you wanted to submit a patch someone would take a  
 look, but the question is whether there is really any calculation  
 whose execution time would be materially affected by this.  Most  
 computations are not dominated by integer power calculations, so is  
 this really worth the trouble?

 Duncan Murdoch



well, the Gnu Scientific Library has the pow_int() functions, which  
are a generalization
of f3(), so someone thinks so. I did a speed test of them but they  
were much slower than
R (for any of f1(), f2(), f3()):

library(gsl)
system.time(ignore - pow_int(a,8))
[1] 1.07 1.11 3.08 0.00 0.00

why the slow execution time?

But I'm far more interested in the philosophy behind your comments.  I
would say that it  definitely *is* worth the trouble because someone,
somewhere, will want fast integer powers, and possibly use R for  
nothing else.

Ken's point about matrix exponentiation is relevant here too.

This is a stated design consideration in Mathematica, I think.

(Of course, I'm not suggesting that other programming tasks be  
suspended!  All I'm pointing
out is that there may exist a user to whom fast integer powers are  
very very important)


very best wishes

rksh


--
Robin Hankin
Uncertainty Analyst
National Oceanography Centre, Southampton
European Way, Southampton SO14 3ZH, UK
  tel  023-8059-7743

__
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] x*x*x*... vs x^n

2005-06-29 Thread Prof Brian Ripley
On Wed, 29 Jun 2005, Duncan Murdoch wrote:

 On 6/29/2005 7:32 AM, Robin Hankin wrote:

 I have been wondering if there one can speed up calculating small powers
 of numbers such as x^8 using multiplication.

 In addition, one can be a bit clever and calculate x^8 using only 3
 multiplies.

 look at this:


  f1 - function(x){x*x*x*x*x*x*x*x}
  f2 - function(x){x^8}
  f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}

 [so f1() and f2() and f3() are algebraically identical]


  a - rnorm(100)
  system.time(ignore - f1(a))
 [1] 0.50 0.17 2.88 0.00 0.00

  system.time(ignore - f2(a))
 [1] 0.31 0.03 1.40 0.00 0.00

  system.time(ignore - f3(a))
 [1] 0.10 0.07 0.18 0.00 0.00


 [these figures show little variance from trial to trial]


 I was expecting f2() and f3() to be about the same.
 I was not expecting a factor of 3 there!

 anyone got any comments?

 If you look in src/main/arithmetic.c, you'll see that R does a general
 real-valued power (using C's pow() function) whenever either one of the
 args is real (except for a few special cases, e.g. non-numbers, or
 powers of 2 or 0.5).  There is an internal R function equivalent to your
 f3, but it is not used in the situation of real^integer (and in any
 case, x^8 is real^real).

 I suppose if you wanted to submit a patch someone would take a look, but
 the question is whether there is really any calculation whose execution
 time would be materially affected by this.  Most computations are not
 dominated by integer power calculations, so is this really worth the
 trouble?

As Luke Tierney frequently points out, selecting special cases can take 
more time than you save.  The assembler code used by modern OSes will
have worked out that compromise pretty effectively: even for real^real it 
spots simple cases of the exponent.

Also, it depends on the CPU: I get on Athlon 2600

 system.time(ignore - f1(a), gcFirst=T)
[1] 0.08 0.05 0.14 0.00 0.00
 system.time(ignore - f2(a), gcFirst=T)
[1] 0.20 0.01 0.20 0.00 0.00
 system.time(ignore - f3(a), gcFirst=T)
[1] 0.03 0.02 0.05 0.00 0.00

and Opteron 248

 system.time(ignore - f1(a), gcFirst=T)
[1] 0.08 0.06 0.14 0.00 0.00
 system.time(ignore - f2(a), gcFirst=T)
[1] 0.19 0.01 0.20 0.00 0.00
 system.time(ignore - f3(a), gcFirst=T)
[1] 0.04 0.01 0.05 0.00 0.00

Note

1) the use of gcFirst=T
2) these need to be run several times to tune the gc() behaviour.  After 
which f1(a) caused a couple of gc()s and the others none.
3) the Opteron is in general a much faster machine so these are all 
surprisingly slow.

Occasionally the C-level pow() is slow: that happened in one MinGW update
and for R Windows users we replaced it.  Even there I was unsure that it 
would make enough difference on a real problem, but eventually found one 
where it did (MDS with Minkowski distances).  It was because of that real 
problem that I added the special case for 0.5, after careful timing
(and hearing Luke's comment alluded to above).

-- 
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] x*x*x*... vs x^n

2005-06-29 Thread Duncan Murdoch
On 6/29/2005 9:31 AM, Robin Hankin wrote:
 Hi  Duncan
 
 
 On Jun 29, 2005, at 02:04 pm, Duncan Murdoch wrote:
 
 On 6/29/2005 7:32 AM, Robin Hankin wrote:

 Hi
 I have been wondering if there one can speed up calculating small  
 powers
 of numbers such as x^8 using multiplication.
 In addition, one can be a bit clever and calculate x^8 using only  
 3  multiplies.
 look at this:
   f1 - function(x){x*x*x*x*x*x*x*x}
   f2 - function(x){x^8}
   f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}

 
 [snip]
 

 If you look in src/main/arithmetic.c, you'll see that R does a  
 general real-valued power (using C's pow() function) whenever  
 either one of the args is real (except for a few special cases,  
 e.g. non-numbers, or powers of 2 or 0.5).  There is an internal R  
 function equivalent to your f3, but it is not used in the situation  
 of real^integer (and in any case, x^8 is real^real).

 I suppose if you wanted to submit a patch someone would take a  
 look, but the question is whether there is really any calculation  
 whose execution time would be materially affected by this.  Most  
 computations are not dominated by integer power calculations, so is  
 this really worth the trouble?

 Duncan Murdoch

 
 
 well, the Gnu Scientific Library has the pow_int() functions, which  
 are a generalization
 of f3(), so someone thinks so. I did a speed test of them but they  
 were much slower than
 R (for any of f1(), f2(), f3()):
 
 library(gsl)
 system.time(ignore - pow_int(a,8))
 [1] 1.07 1.11 3.08 0.00 0.00
 
 why the slow execution time?

Shouldn't you ask the gsl maintainer that? :-)

 But I'm far more interested in the philosophy behind your comments.  I
 would say that it  definitely *is* worth the trouble because someone,
 somewhere, will want fast integer powers, and possibly use R for  
 nothing else.
 
 Ken's point about matrix exponentiation is relevant here too.
 
 This is a stated design consideration in Mathematica, I think.
 
 (Of course, I'm not suggesting that other programming tasks be  
 suspended!  All I'm pointing
 out is that there may exist a user to whom fast integer powers are  
 very very important)

Then that user should submit the patch, not you.  But whoever does it 
should include an argument to convince an R core member that the change 
is worth looking at, and do it well enough that the patch is accepted.
Changes to the way R does arithmetic affect everyone, so they had better 
be done right, and checking them takes time.

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] all connections are in use error during lazyload stage of packa ge installation

2005-06-29 Thread Tuszynski, Jaroslaw W.
Hi,

I suddenly started getting strange errors while working on my caTools
package:

 RCMD install C:/programs/R/rw2011/src/library/caTools
  ..
preparing package caTools for lazy loading
Error in file(file, r, encoding = encoding) :
all connections are in use
Execution halted
make: *** [lazyload] Error 1
*** Installation of caTools failed ***

I searched through R website and mailing lists but did not found many clues.
I am not running R environment so showConnections and CloseAllConnections
functions are not helpful. Did I run into some OS specific limit on number
of files allowed? I closed all other programs and rebooted my machine but it
did not help.

Any hints would be appreciated.

My System is:
- R version: R 2.1.1
- Operating System: Win XP
- Compiler: mingw32-gcc-3.4.2

Jarek
\===

 Jarek Tuszynski, PhD.   o / \ 
 Science Applications International Corporation  \__,|  
 (703) 676-4192  \
 [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] all connections are in use error during lazyload stage of packa ge installation

2005-06-29 Thread Gabor Grothendieck
On 6/29/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote:
 Hi,
 
 I suddenly started getting strange errors while working on my caTools
 package:
 
 RCMD install C:/programs/R/rw2011/src/library/caTools
  ..
preparing package caTools for lazy loading
Error in file(file, r, encoding = encoding) :
all connections are in use
Execution halted
make: *** [lazyload] Error 1
*** Installation of caTools failed ***
 
 I searched through R website and mailing lists but did not found many clues.
 I am not running R environment so showConnections and CloseAllConnections
 functions are not helpful. Did I run into some OS specific limit on number
 of files allowed? I closed all other programs and rebooted my machine but it
 did not help.

I don't know what the problem is but you could see if adding

LazyLoad: no

to the DESCRIPTION file causes it to go away.  Even if it does it
would probably be best to track down what the cause is and in the past 
when I have had problems running R CMD on a package of mine I 
repeatedly removed half my package and reran R CMD until I found the 
offending portion.

__
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] x*x*x*... vs x^n

2005-06-29 Thread Robin Hankin

On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote:

 On 6/29/2005 9:31 AM, Robin Hankin wrote:

 Hi  Duncan


 library(gsl)
 system.time(ignore - pow_int(a,8))
 [1] 1.07 1.11 3.08 0.00 0.00

 why the slow execution time?


 Shouldn't you ask the gsl maintainer that? :-)


well  I did ask myself, but  in this case the gsl maintainer
told me to ask the gsl co-author, who
is no doubt much better informed in these matters ;-)


 (Of course, I'm not suggesting that other programming tasks be
 suspended!  All I'm pointing
 out is that there may exist a user to whom fast integer powers are
 very very important)


 Then that user should submit the patch, not you.  But whoever does it
 should include an argument to convince an R core member that the  
 change
 is worth looking at, and do it well enough that the patch is accepted.
 Changes to the way R does arithmetic affect everyone, so they had  
 better
 be done right, and checking them takes time.


yes, that's a fair point.
But including a native R command pow.int(), say, wouldn't  affect  
anyone, would it?
One could even use the (tested) GSL code, as it is GPL'ed.

This would just be a new function that users could use at their  
discretion, and no
existing code would break.

I assume that such a function would  not suffer whatever performance  
disadvantage
that the GSL package approach had, so it may well be quite a  
significant improvement
over the method used by R_pow_di() in arithmetic.c especially for  
large n.


best wishes

rksh

 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-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] moving correlation coef ?

2005-06-29 Thread vincent
Hello,

R gives us the correlation functions cor(). (Many thanks ;-))
Does it also exist a moving correlation coefficient ?
(like the moving average).
If not, could someone give me some infos or link
on how to practically implement such a function in R.

(I did search for moving correlation on the R homepage
but didn't find anything.)

Thank you.
Vincent

__
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] moving correlation coef ?

2005-06-29 Thread Sundar Dorai-Raj


vincent wrote:
 Hello,
 
 R gives us the correlation functions cor(). (Many thanks ;-))
 Does it also exist a moving correlation coefficient ?
 (like the moving average).
 If not, could someone give me some infos or link
 on how to practically implement such a function in R.
 
 (I did search for moving correlation on the R homepage
 but didn't find anything.)
 
 Thank you.
 Vincent
 


See ?running in the gtools package:

library(gtools)
X - rnorm(100); Y - rnorm(100)
running(X, Y, cor)

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] all connections are in use error during lazyload stage of packa ge installation

2005-06-29 Thread Tuszynski, Jaroslaw W.
 I found the problem, by doing comparison of directories and files of
working and not working versions, and applying changes one by one until one
caused install to fail. It was a case of having a call to a source
function somewhere in my code, that I forgot about. 

I was definitely doing something silly , but it was not a meaningful error
message.

Jarek
\===

 Jarek Tuszynski, PhD.   o / \ 
 Science Applications International Corporation  \__,|  
 (703) 676-4192  \
 [EMAIL PROTECTED] `\



-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] 
Sent: Wednesday, June 29, 2005 10:07 AM
To: Tuszynski, Jaroslaw W.
Cc: r help
Subject: Re: [R] all connections are in use error during lazyload stage of
packa ge installation

On 6/29/05, Tuszynski, Jaroslaw W. [EMAIL PROTECTED] wrote:
 Hi,
 
 I suddenly started getting strange errors while working on my caTools
 package:
 
 RCMD install C:/programs/R/rw2011/src/library/caTools
  ..
preparing package caTools for lazy loading
Error in file(file, r, encoding = encoding) :
all connections are in use
Execution halted
make: *** [lazyload] Error 1
*** Installation of caTools failed ***
 
 I searched through R website and mailing lists but did not found many
clues.
 I am not running R environment so showConnections and 
 CloseAllConnections functions are not helpful. Did I run into some OS 
 specific limit on number of files allowed? I closed all other programs 
 and rebooted my machine but it did not help.

I don't know what the problem is but you could see if adding

LazyLoad: no

to the DESCRIPTION file causes it to go away.  Even if it does it would
probably be best to track down what the cause is and in the past when I have
had problems running R CMD on a package of mine I repeatedly removed half my
package and reran R CMD until I found the offending portion.

__
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] x*x*x*... vs x^n

2005-06-29 Thread Ravi Varadhan

I ran 100 repetitions of the 3 multiplications that Robin had compared.
Here are the summaries of system times (I only took the first component of
system.time) that I obtained.  It is clear that f1() is nearly twice as slow
as f2() which is slightly slower (not 3 times slower as claimed by Robin)
than f3(). So, I don't think that there is much to choose between the
cleverer way and the most obvious way to compute integer powers.

 summary(f1time)
   Min. 1st Qu.  MedianMean 3rd Qu.Max. 
  0.060   0.170   0.210   0.199   0.230   0.300 
 summary(f2time)
   Min. 1st Qu.  MedianMean 3rd Qu.Max. 
  0.060   0.100   0.110   0.128   0.170   0.190 
 summary(f3time)
   Min. 1st Qu.  MedianMean 3rd Qu.Max. 
 0.  0.0300  0.0950  0.0779  0.1100  0.1300 

Ravi.

--
Ravi Varadhan, Ph.D.
Assistant Professor,  The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email:  [EMAIL PROTECTED]
--
 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help-
 [EMAIL PROTECTED] On Behalf Of Tuszynski, Jaroslaw W.
 Sent: Wednesday, June 29, 2005 8:25 AM
 To: 'Robin Hankin'; r-help
 Subject: Re: [R] x*x*x*... vs x^n
 
 I tried your code and got different results:
system.time(ignore - f1(a))
   [1] 0.83 0.09 1.08   NA   NA
system.time(ignore - f2(a))
   [1] 0.38 0.01 0.41   NA   NA
system.time(ignore - f3(a))
   [1] 0.32 0.04 0.43   NA   NA
 
 So I tried it again but with a loop and got:
 
 for(i in 1:10) cat(system.time(ignore - f2(a)), \n)
   0.36 0.04 0.44 NA NA
   0.32 0.01 0.34 NA NA
   0.28 0.03 0.32 NA NA
   0.29 0.03 0.35 NA NA
   0.3 0.02 0.32 NA NA
   0.28 0.03 0.32 NA NA
   0.3 0.02 0.32 NA NA
   0.29 0.02 0.34 NA NA
   0.23 0.03 0.32 NA NA
   0.42 0 0.45 NA NA
 
for(i in 1:10) cat(system.time(ignore - f3(a)), \n)
   0.19 0.04 0.25 NA NA
   0.17 0.04 0.25 NA NA
   0.21 0.02 0.25 NA NA
   0.21 0.02 0.23 NA NA
   0.18 0.04 0.23 NA NA
   0.18 0.05 0.23 NA NA
   0.18 0.04 0.25 NA NA
   0.17 0.06 0.23 NA NA
   0.2 0.02 0.23 NA NA
   0.14 0.06 0.25 NA NA
 
 It seems to me that f3 is 50% slower than f2 not 300%.
 
 My System is:
 - R version: R 2.1.1
 - Operating System: Win XP
 - Compiler: mingw32-gcc-3.4.2
 
 Jarek
 \===
 
  Jarek Tuszynski, PhD.   o / \
  Science Applications International Corporation  \__,|
  (703) 676-4192  \
  [EMAIL PROTECTED] `\
 
 
 
 -Original Message-
 From: [EMAIL PROTECTED]
 [mailto:[EMAIL PROTECTED] On Behalf Of Robin Hankin
 Sent: Wednesday, June 29, 2005 7:32 AM
 To: r-help
 Subject: [R] x*x*x*... vs x^n
 
 Hi
 
 I have been wondering if there one can speed up calculating small powers
 of
 numbers such as x^8 using multiplication.
 
 In addition, one can be a bit clever and calculate x^8 using only 3
 multiplies.
 
 look at this:
 
 
   f1 - function(x){x*x*x*x*x*x*x*x}
   f2 - function(x){x^8}
   f3 - function(x){x2 - x*x;x4 - x2*x2;return(x4*x4)}
 
 [so f1() and f2() and f3() are algebraically identical]
 
 
   a - rnorm(100)
   system.time(ignore - f1(a))
 [1] 0.50 0.17 2.88 0.00 0.00
 
   system.time(ignore - f2(a))
 [1] 0.31 0.03 1.40 0.00 0.00
 
   system.time(ignore - f3(a))
 [1] 0.10 0.07 0.18 0.00 0.00
 
 
 [these figures show little variance from trial to trial]
 
 
 I was expecting f2() and f3() to be about the same.
 I was not expecting a factor of 3 there!
 
 anyone got any comments?
 
 
 
 
 --
 Robin Hankin
 Uncertainty Analyst
 National Oceanography Centre, Southampton European Way, Southampton SO14
 3ZH, UK
   tel  023-8059-7743
 
 __
 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-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] x*x*x*... vs x^n

2005-06-29 Thread davidr
In general, the Russian peasant algorithm, which requires only O(log
n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of
Computer Programming. Volume 2: Seminumerical Algorithms has an in depth
discussion.
I have had to use this in the past, when computers were slower and
compilers were not so clever. It is also better when x is not just a
real number, say complex or matrix (as has been mentioned.)
In many cases though, one needs many powers sequentially, and then it
may be more efficient to just multiply the previous power by x and use
the power, etc. (unless you have a parallel computer.)
So 

pows - x^(1:1000) 
# use pows in calculations

could be sped up by employing a faster algorithm, but probably a loop
will be faster:

pows - 1
for(i in 1:1000) {
  pows - pows * x
  # use this power
}

David L. Reiner, Ph.D.
Rho Trading
 

 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help-
 [EMAIL PROTECTED] On Behalf Of Robin Hankin
 Sent: Wednesday, June 29, 2005 9:13 AM
 To: Duncan Murdoch
 Cc: r-help; Robin Hankin
 Subject: Re: [R] x*x*x*... vs x^n
 
 
 On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote:
 
  On 6/29/2005 9:31 AM, Robin Hankin wrote:
 
  Hi  Duncan
 
 
  library(gsl)
  system.time(ignore - pow_int(a,8))
  [1] 1.07 1.11 3.08 0.00 0.00
 
  why the slow execution time?
 
 
  Shouldn't you ask the gsl maintainer that? :-)
 
 
 well  I did ask myself, but  in this case the gsl maintainer
 told me to ask the gsl co-author, who
 is no doubt much better informed in these matters ;-)
 
 
  (Of course, I'm not suggesting that other programming tasks be
  suspended!  All I'm pointing
  out is that there may exist a user to whom fast integer powers are
  very very important)
 
 
  Then that user should submit the patch, not you.  But whoever does
it
  should include an argument to convince an R core member that the
  change
  is worth looking at, and do it well enough that the patch is
accepted.
  Changes to the way R does arithmetic affect everyone, so they had
  better
  be done right, and checking them takes time.
 
 
 yes, that's a fair point.
 But including a native R command pow.int(), say, wouldn't  affect
 anyone, would it?
 One could even use the (tested) GSL code, as it is GPL'ed.
 
 This would just be a new function that users could use at their
 discretion, and no
 existing code would break.
 
 I assume that such a function would  not suffer whatever performance
 disadvantage
 that the GSL package approach had, so it may well be quite a
 significant improvement
 over the method used by R_pow_di() in arithmetic.c especially for
 large n.
 
 
 best wishes
 
 rksh
 
  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-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] Selecting rows regarding the frequency of a factor variable.

2005-06-29 Thread Ghislain Vieilledent
Hi and sorry to disturb,

I'll try to be as clear as possible:
I want to select rows of a data frame called Data2.Iso regarding the 
frequency of a factor variable called Variete that I want =4.

I used function table to have the frequency:
  FRAMEVARIETE-as.data.frame(table(Data2.Iso$Variete))
Then I selected the modalities with a frequency =4:
  FRAMEVARIETE2-FRAMEVARIETE[FRAMEVARIETE$Freq=4,]
  as.character(FRAMEVARIETE2[,Variete])

But then, how to select the rows with those modalities?
Does anyone can help me?

Thanks!

Ghislain.

__
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] moving correlation coef ?

2005-06-29 Thread Spencer Graves
or ?rapply in package zoo.

  spencer graves

Sundar Dorai-Raj wrote:

 
 vincent wrote:
 
Hello,

R gives us the correlation functions cor(). (Many thanks ;-))
Does it also exist a moving correlation coefficient ?
(like the moving average).
If not, could someone give me some infos or link
on how to practically implement such a function in R.

(I did search for moving correlation on the R homepage
but didn't find anything.)

Thank you.
Vincent

 
 
 
 See ?running in the gtools package:
 
 library(gtools)
 X - rnorm(100); Y - rnorm(100)
 running(X, Y, cor)
 
 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

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


Re: [R] Selecting rows regarding the frequency of a factor variab le.

2005-06-29 Thread Liaw, Andy
See if this does what you want:

 dat - data.frame(f=factor(sample(letters[1:10], 100, replace=TRUE)),
x=runif(100))
 str(dat)
`data.frame':   100 obs. of  2 variables:
 $ f: Factor w/ 10 levels a,b,c,d,..: 2 5 10 9 10 3 9 8 3 1 ...
 $ x: num  0.9162 0.0481 0.3048 0.0938 0.8599 ...
 g - names(which(table(dat$f)  11))
 g
[1] c j
 dat[dat$f %in% g,]
   f  x
3  j 0.30477413
5  j 0.85992597
6  c 0.86881528
9  c 0.87317095
16 c 0.84252048
18 j 0.24039606
19 j 0.58927414
21 j 0.10077745
32 j 0.72275870
35 c 0.26001549
37 j 0.09608521
40 c 0.15481625
44 c 0.70203309
47 c 0.95907223
50 j 0.35258966
54 c 0.93422614
58 c 0.36546841
61 c 0.55123183
64 j 0.82995122
65 c 0.89104229
66 j 0.81661377
77 j 0.21134708
87 c 0.16602335
92 c 0.02175573
96 j 0.97864088

Andy

 From:  Ghislain Vieilledent
 
 
 Hi and sorry to disturb,
 
 I'll try to be as clear as possible:
 I want to select rows of a data frame called Data2.Iso 
 regarding the 
 frequency of a factor variable called Variete that I want =4.
 
 I used function table to have the frequency:
   FRAMEVARIETE-as.data.frame(table(Data2.Iso$Variete))
 Then I selected the modalities with a frequency =4:
   FRAMEVARIETE2-FRAMEVARIETE[FRAMEVARIETE$Freq=4,]
   as.character(FRAMEVARIETE2[,Variete])
 
 But then, how to select the rows with those modalities?
 Does anyone can help me?
 
 Thanks!
 
 Ghislain.
 
 __
 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] MLE with optim

2005-06-29 Thread Carsten Steinhoff
Hello,
 
I tried to fit a lognormal distribution by using optim. But sadly the output
seems to be incorrect.
Who can tell me where the bug is?
 
test = rlnorm(100,5,3)
logL= function(parm, x,...) -sum(log(dlnorm(x,parm,...)))
start= list(meanlog=5, sdlog=3)
optim(start,logL,x=test)$par
 
Carsten.

[[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] x*x*x*... vs x^n

2005-06-29 Thread davidr
Looking at the code for gsl_pow_int, I see they do use that method.

David L. Reiner
 

 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help-
 [EMAIL PROTECTED] On Behalf Of David Reiner
 [EMAIL PROTECTED]
 Sent: Wednesday, June 29, 2005 9:50 AM
 To: r-help
 Subject: Re: [R] x*x*x*... vs x^n
 
 In general, the Russian peasant algorithm, which requires only O(log
 n) multiplications, is very good. Section 4.6.3 of Knuth's The Art of
 Computer Programming. Volume 2: Seminumerical Algorithms has an in
depth
 discussion.
 I have had to use this in the past, when computers were slower and
 compilers were not so clever. It is also better when x is not just a
 real number, say complex or matrix (as has been mentioned.)
 In many cases though, one needs many powers sequentially, and then it
 may be more efficient to just multiply the previous power by x and use
 the power, etc. (unless you have a parallel computer.)
 So
 
 pows - x^(1:1000)
 # use pows in calculations
 
 could be sped up by employing a faster algorithm, but probably a loop
 will be faster:
 
 pows - 1
 for(i in 1:1000) {
   pows - pows * x
   # use this power
 }
 
 David L. Reiner, Ph.D.
 Rho Trading
 
 
  -Original Message-
  From: [EMAIL PROTECTED] [mailto:r-help-
  [EMAIL PROTECTED] On Behalf Of Robin Hankin
  Sent: Wednesday, June 29, 2005 9:13 AM
  To: Duncan Murdoch
  Cc: r-help; Robin Hankin
  Subject: Re: [R] x*x*x*... vs x^n
 
 
  On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote:
 
   On 6/29/2005 9:31 AM, Robin Hankin wrote:
  
   Hi  Duncan
  
  
   library(gsl)
   system.time(ignore - pow_int(a,8))
   [1] 1.07 1.11 3.08 0.00 0.00
  
   why the slow execution time?
  
  
   Shouldn't you ask the gsl maintainer that? :-)
  
 
  well  I did ask myself, but  in this case the gsl maintainer
  told me to ask the gsl co-author, who
  is no doubt much better informed in these matters ;-)
 
  
   (Of course, I'm not suggesting that other programming tasks be
   suspended!  All I'm pointing
   out is that there may exist a user to whom fast integer powers
are
   very very important)
  
  
   Then that user should submit the patch, not you.  But whoever does
 it
   should include an argument to convince an R core member that the
   change
   is worth looking at, and do it well enough that the patch is
 accepted.
   Changes to the way R does arithmetic affect everyone, so they had
   better
   be done right, and checking them takes time.
  
 
  yes, that's a fair point.
  But including a native R command pow.int(), say, wouldn't  affect
  anyone, would it?
  One could even use the (tested) GSL code, as it is GPL'ed.
 
  This would just be a new function that users could use at their
  discretion, and no
  existing code would break.
 
  I assume that such a function would  not suffer whatever performance
  disadvantage
  that the GSL package approach had, so it may well be quite a
  significant improvement
  over the method used by R_pow_di() in arithmetic.c especially for
  large n.
 
 
  best wishes
 
  rksh
 
   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-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-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] MLE with optim

2005-06-29 Thread Sundar Dorai-Raj


Carsten Steinhoff wrote:
 Hello,
  
 I tried to fit a lognormal distribution by using optim. But sadly the output
 seems to be incorrect.
 Who can tell me where the bug is?
  
 test = rlnorm(100,5,3)
 logL= function(parm, x,...) -sum(log(dlnorm(x,parm,...)))
 start= list(meanlog=5, sdlog=3)
 optim(start,logL,x=test)$par
  
 Carsten.
 

You are only supplying the meanlog argument to dlnorm. Also you can set 
log = TRUE in dlnorm to avoid the log computation after calling dlnorm.

set.seed(1)
x - rlnorm(100, 5, 3)
logL - function(par, x) -sum(dlnorm(x, par[1], par[2], TRUE))
start - list(meanlog = 5, sdlog = 2)
optim(start, logL, x = x)

Why not just use MASS::fitdistr instead?

library(MASS)
fitdistr(x, dlnorm, start)

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] MLE with optim

2005-06-29 Thread Dimitris Rizopoulos
the following work for me:


x - rlnorm(1000, 5, 3)

fn - function(parms, dat) -sum(dlnorm(dat, parms[1], parms[2], log = 
TRUE))
optim(c(5, 3), fn, dat = x)

library(MASS)
fitdistr(x, log-normal, list(meanlog = 5,  sdlog = 3))


I hope it helps.

Best,
Dimitris


Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.be/biostat/
 http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm


- Original Message - 
From: Carsten Steinhoff [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Sent: Wednesday, June 29, 2005 5:19 PM
Subject: [R] MLE with optim


 Hello,

 I tried to fit a lognormal distribution by using optim. But sadly 
 the output
 seems to be incorrect.
 Who can tell me where the bug is?

 test = rlnorm(100,5,3)
 logL= function(parm, x,...) -sum(log(dlnorm(x,parm,...)))
 start= list(meanlog=5, sdlog=3)
 optim(start,logL,x=test)$par

 Carsten.

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


Re: [R] moving correlation coef ?

2005-06-29 Thread vincent
Thank you for your answers.

In fact, i believe my question wasn't precise enough.
I don't want to have a moving/sliding windows over the data
to correlate (i am already doing that).

If I have 2 vectors
X = (x1, x2, x3, ..., xt)
Y = (y1, y2, x3, ..., yt)
I want the most recent elements (t) to have a heavier weight
in the correlation calculus than the older ones (1).

I think that one simple way to do that is for example to
compute cor() over XP, YP where
XP = (x1, x2, x2, x3, x3, x3, ..., xt, xt, xt)
YP = (y1, y2, y2, y3, y3, y3, ..., yt, yt, yt)
ie where each element is repeated several times according
to its freshness.
It's quite naive ! so if there is a cleaver idea, many thanks.

Thanks
Vincent

__
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] poly() in lm() leads to wrong coefficients (but correct residuals)

2005-06-29 Thread Andreas Neumann
Dear all,

I am using poly() in lm() in the following form.

1 DelsDPWOS.lm3 - lm(DelsPDWOS[,1] ~ poly(DelsPDWOS[,4],3))

2 DelsDPWOS.I.lm3 - lm(DelsPDWOS[,1] ~ poly(I(DelsPDWOS[,4]),3))

3 DelsDPWOS.2.lm3 -
   lm(DelsPDWOS[,1]~DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3))

1 and 2 lead to identical but wrong results. 3 is correct. Surprisingly
(to me) the residuals are the same (correct) in all cases although the
coefficients of 1 (and 2) are wrong and do not in any way match the
stated regression polynomial. (summaries below)

QUESTION:
Is there a correct way to use poly() in lm() since version 3 becomes quite
tedious if used more often especially with higher order polynomials?

I am using R Version 1.9.0 (2004-04-12).

Thank you for your time,

Andreas Neumann



SUMMARIES:
 summary(DelsDPWOS.lm3)

Call:
lm(formula = DelsPDWOS[, 1] ~ poly(DelsPDWOS[, 4], 3))

Residuals:
   12345678
-0.11414  0.12756 -0.21060  0.04636 -0.03244  0.16030  0.04290 -0.03944
   9
 0.01949

Coefficients:
 Estimate Std. Error t value Pr(|t|)
(Intercept)   1.30.04861  27.430 1.21e-06 ***
poly(DelsPDWOS[, 4], 3)1 -1.274640.14583  -8.741 0.000325 ***
poly(DelsPDWOS[, 4], 3)2  0.274830.14583   1.885 0.118175
poly(DelsPDWOS[, 4], 3)3  0.115900.14583   0.795 0.462774
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 0.1458 on 5 degrees of freedom
Multiple R-Squared: 0.9416, Adjusted R-squared: 0.9065
F-statistic: 26.86 on 3 and 5 DF,  p-value: 0.001645

 summary(DelsDPWOS.2.lm3)

Call:
lm(formula = DelsPDWOS[, 1] ~ DelsPDWOS[, 4] + I(DelsPDWOS[,
4]^2) + I(DelsPDWOS[, 4]^3))

Residuals:
   12345678
-0.11414  0.12756 -0.21060  0.04636 -0.03244  0.16030  0.04290 -0.03944
   9
 0.01949

Coefficients:
Estimate Std. Error t value Pr(|t|)
(Intercept)   -1.486  7.242  -0.2050.846
DelsPDWOS[, 4] 8.640 14.161   0.6100.568
I(DelsPDWOS[, 4]^2)   -6.526  8.799  -0.7420.492
I(DelsPDWOS[, 4]^3)1.375  1.730   0.7950.463

Residual standard error: 0.1458 on 5 degrees of freedom
Multiple R-Squared: 0.9416, Adjusted R-squared: 0.9065
F-statistic: 26.86 on 3 and 5 DF,  p-value: 0.001645

__
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] poly() in lm() leads to wrong coefficients (but correct residuals)

2005-06-29 Thread Markus Jäntti
On Wed, 2005-06-29 at 18:19 +0200, Andreas Neumann wrote:
 Dear all,
 
 I am using poly() in lm() in the following form.
 
 1 DelsDPWOS.lm3 - lm(DelsPDWOS[,1] ~ poly(DelsPDWOS[,4],3))
 
 2 DelsDPWOS.I.lm3 - lm(DelsPDWOS[,1] ~ poly(I(DelsPDWOS[,4]),3))
 
 3 DelsDPWOS.2.lm3 -
lm(DelsPDWOS[,1]~DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3))
 
 1 and 2 lead to identical but wrong results. 3 is correct. Surprisingly
 (to me) the residuals are the same (correct) in all cases although the
 coefficients of 1 (and 2) are wrong and do not in any way match the
 stated regression polynomial. (summaries below)
 
 QUESTION:
 Is there a correct way to use poly() in lm() since version 3 becomes quite
 tedious if used more often especially with higher order polynomials?
 

The coefficients  using 1 and 2 are not incorrect.  
poly() defines orthogonal polynomials, whereas your
DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3 
contruct defines an ordinary polynomial. 

You should be able to recover version 3 coefficients from 1 and 2.
See help(poly)

 x - runif(10)
 x
 [1] 0.1878 0.2415 0.5834 0.6556 0.4112 0.3399 0.8144 0.1134 0.7360
0.0463
 model.matrix(~ poly(x, 2))
   (Intercept) poly(x, 2)1 poly(x, 2)2
11-0.27648 -0.0452
21-0.21052 -0.1899
31 0.20937 -0.2708
41 0.29799 -0.1021
51-0.00212 -0.4117
61-0.08970 -0.3621
71 0.49297  0.4968
81-0.36790  0.2148
91 0.39672  0.1620
10   1-0.45033  0.5082
attr(,assign)
[1] 0 1 1
 model.matrix(~ x + I(x^2))
   (Intercept)  x  I(x^2)
11 0.1878 0.03528
21 0.2415 0.05834
31 0.5834 0.34040
41 0.6556 0.42982
51 0.4112 0.16911
61 0.3399 0.11554
71 0.8144 0.66320
81 0.1134 0.01286
91 0.7360 0.54169
10   1 0.0463 0.00214
attr(,assign)
[1] 0 1 2



Regards,

-- 
Markus Jantti
Abo Akademi University
[EMAIL PROTECTED]
http://www.iki.fi/~mjantti

__
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] moving correlation coef ?

2005-06-29 Thread Sundar Dorai-Raj


vincent wrote:
 Thank you for your answers.
 
 In fact, i believe my question wasn't precise enough.
 I don't want to have a moving/sliding windows over the data
 to correlate (i am already doing that).
 
 If I have 2 vectors
 X = (x1, x2, x3, ..., xt)
 Y = (y1, y2, x3, ..., yt)
 I want the most recent elements (t) to have a heavier weight
 in the correlation calculus than the older ones (1).
 
 I think that one simple way to do that is for example to
 compute cor() over XP, YP where
 XP = (x1, x2, x2, x3, x3, x3, ..., xt, xt, xt)
 YP = (y1, y2, y2, y3, y3, y3, ..., yt, yt, yt)
 ie where each element is repeated several times according
 to its freshness.
 It's quite naive ! so if there is a cleaver idea, many thanks.
 
 Thanks
 Vincent
 

Perhaps ?cov.wt will work for you? Your example would be identical to:

set.seed(1)
X - rnorm(100); Y - rnorm(100)
# using cov.wt
rho1 - cov.wt(cbind(X, Y), 1:100, cor = TRUE)$cor[1, 2]
# your weighting scheme
rho2 - cor(X[rep(1:100, 1:100)], Y[rep(1:100, 1:100)])
all.equal(rho1, rho2)
# [1] TRUE

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] moving correlation coef ?

2005-06-29 Thread davidr

One common weighting scheme is exponentially weighted, i.e., wt =
L^(0:m) ,
where 0  L = 1 .

David L. Reiner

p.s.
If your question is coming from a financial application, you might be
interested in the R-sig-finance list, as well as reading the RiskMetrics
(r)
document Return to RiskMetrics: The Evolution of a Standard, by Jorge
Mina and Jerry Yi Xiao (available at their site, after a free
registration) where exponentially weighted moving statistics are
discussed at length.

 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help-
 [EMAIL PROTECTED] On Behalf Of vincent
 Sent: Wednesday, June 29, 2005 11:02 AM
 Cc: [EMAIL PROTECTED]
 Subject: Re: [R] moving correlation coef ?
 
 Thank you for your answers.
 
 In fact, i believe my question wasn't precise enough.
 I don't want to have a moving/sliding windows over the data
 to correlate (i am already doing that).
 
 If I have 2 vectors
 X = (x1, x2, x3, ..., xt)
 Y = (y1, y2, x3, ..., yt)
 I want the most recent elements (t) to have a heavier weight
 in the correlation calculus than the older ones (1).
 
 I think that one simple way to do that is for example to
 compute cor() over XP, YP where
 XP = (x1, x2, x2, x3, x3, x3, ..., xt, xt, xt)
 YP = (y1, y2, y2, y3, y3, y3, ..., yt, yt, yt)
 ie where each element is repeated several times according
 to its freshness.
 It's quite naive ! so if there is a cleaver idea, many thanks.
 
 Thanks
 Vincent
 
 __
 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] Help with regression modeling

2005-06-29 Thread Tony Young
Hello all,

 

I'm using R version 2.0.1. I have been having trouble with my linear
modeling. I have a table that looks something like this:

 

T  RSS DSS  LSPFCOLS   PS
R  RTT  Actual   Max  COMM

char  5   MSS 2  1251  1
4096 450   0.00164.300.13  64.17

char  5   MSS 50 1251  25
4096 450   0.00174.970.14  74.83

char  5   MSS 200   1251  100
4096 450   0.001111.12  0.13  110.99

char  5   MSS 3  1251  1
4096 450   0.00164.310.14  64.17

char  5   MSS 75 1251  25
4096 450   0.00175.990.14  75.85

char  5   MSS 300   1251  100
4096 450   0.001116.62  0.14  116.48

 

that I was modeling using the R command

 

model.f - lm(comm ~ I(rtt * ceiling(rss/pf)) + I((s * rss) / ls) +
ds*(I(((s+s^2)*rss)/r) + I(((s+s^2)*rss)/(r*ps

 

R gives me a table like this for coefficients

 

  I(rtt * ceiling(rss/pf))I((s *
rss)/ls) 

  8.174598e-01
1.353734e+00 

 dsASA
dsDB2 

  3.101663e+00
2.587372e+00 

 dsMSS
dsOracle 

  1.575822e+01
9.092041e+00 

I(((s + s^2) * rss)/r)  I(((s + s^2) * rss)/(r *
ps)) 

  2.151170e-07
3.087502e-05 

  dsDB2:I(((s + s^2) * rss)/r)   dsMSS:I(((s + s^2) *
rss)/r) 

 -6.972409e-08
1.812149e-07 

   dsOracle:I(((s + s^2) * rss)/r)dsDB2:I(((s + s^2) * rss)/(r *
ps)) 

  1.251699e-08
6.687324e-05 

   dsMSS:I(((s + s^2) * rss)/(r * ps)) dsOracle:I(((s + s^2) * rss)/(r *
ps)) 

  4.731049e-05
-9.243794e-05

 

What I want at the end of the day is a formula that says given this set of
parameters, here is the predicted value. When I plug the above coefficients
into Excel, I get a value that is WAY off the fitted value that R gives me.

 

My questions are:

 

1) How do I use these coefficients in a formula?

 

2) How can I construct a formula from these coefficients that I can publish
in my thesis?

 

Thanks for your help!

 

ty


[[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] x*x*x*... vs x^n

2005-06-29 Thread davidr
I was surprised to find that I was wrong about powers of complexes:

 seq.pow1 - function(x,n) {
+   y - rep(x,n)
+   for(i in 2:n) y[i] - y[i-1] * x
+   y
+ }
 seq.pow2 - function(x,n) x^(1:n)
 x - 1.001 + 1i * 0.999
# several reps of the following
 system.time(ignore - seq.pow1(x,10),gcFirst=TRUE)
[1] 0.73 0.00 0.74   NA   NA
# several reps of the following
 system.time(ignore - seq.pow2(x,10),gcFirst=TRUE)
[1] 0.35 0.00 0.35   NA   NA

I apologize for using probably below when not was correct (modulo
grammar).

David L. Reiner
 

 -Original Message-
 From: [EMAIL PROTECTED] [mailto:r-help-
 [EMAIL PROTECTED] On Behalf Of David Reiner
 [EMAIL PROTECTED]
 Sent: Wednesday, June 29, 2005 10:24 AM
 To: r-help
 Subject: Re: [R] x*x*x*... vs x^n
 
 Looking at the code for gsl_pow_int, I see they do use that method.
 
 David L. Reiner
 
 
  -Original Message-
  From: [EMAIL PROTECTED] [mailto:r-help-
  [EMAIL PROTECTED] On Behalf Of David Reiner
  [EMAIL PROTECTED]
  Sent: Wednesday, June 29, 2005 9:50 AM
  To: r-help
  Subject: Re: [R] x*x*x*... vs x^n
 
  In general, the Russian peasant algorithm, which requires only
O(log
  n) multiplications, is very good. Section 4.6.3 of Knuth's The Art
of
  Computer Programming. Volume 2: Seminumerical Algorithms has an in
 depth
  discussion.
  I have had to use this in the past, when computers were slower and
  compilers were not so clever. It is also better when x is not just a
  real number, say complex or matrix (as has been mentioned.)
  In many cases though, one needs many powers sequentially, and then
it
  may be more efficient to just multiply the previous power by x and
use
  the power, etc. (unless you have a parallel computer.)
  So
 
  pows - x^(1:1000)
  # use pows in calculations
 
  could be sped up by employing a faster algorithm, but probably a
loop
  will be faster:
 
  pows - 1
  for(i in 1:1000) {
pows - pows * x
# use this power
  }
 
  David L. Reiner, Ph.D.
  Rho Trading
 
 
   -Original Message-
   From: [EMAIL PROTECTED] [mailto:r-help-
   [EMAIL PROTECTED] On Behalf Of Robin Hankin
   Sent: Wednesday, June 29, 2005 9:13 AM
   To: Duncan Murdoch
   Cc: r-help; Robin Hankin
   Subject: Re: [R] x*x*x*... vs x^n
  
  
   On Jun 29, 2005, at 02:47 pm, Duncan Murdoch wrote:
  
On 6/29/2005 9:31 AM, Robin Hankin wrote:
   
Hi  Duncan
   
   
library(gsl)
system.time(ignore - pow_int(a,8))
[1] 1.07 1.11 3.08 0.00 0.00
   
why the slow execution time?
   
   
Shouldn't you ask the gsl maintainer that? :-)
   
  
   well  I did ask myself, but  in this case the gsl maintainer
   told me to ask the gsl co-author, who
   is no doubt much better informed in these matters ;-)
  
   
(Of course, I'm not suggesting that other programming tasks be
suspended!  All I'm pointing
out is that there may exist a user to whom fast integer powers
 are
very very important)
   
   
Then that user should submit the patch, not you.  But whoever
does
  it
should include an argument to convince an R core member that the
change
is worth looking at, and do it well enough that the patch is
  accepted.
Changes to the way R does arithmetic affect everyone, so they
had
better
be done right, and checking them takes time.
   
  
   yes, that's a fair point.
   But including a native R command pow.int(), say, wouldn't  affect
   anyone, would it?
   One could even use the (tested) GSL code, as it is GPL'ed.
  
   This would just be a new function that users could use at their
   discretion, and no
   existing code would break.
  
   I assume that such a function would  not suffer whatever
performance
   disadvantage
   that the GSL package approach had, so it may well be quite a
   significant improvement
   over the method used by R_pow_di() in arithmetic.c especially for
   large n.
  
  
   best wishes
  
   rksh
  
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-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-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 

[R] Generalized Linear Mixed Models

2005-06-29 Thread Francisco . Redelico
Hello! 

I am trying to fit a Generalized Linear Mixed Model, ordinally I use 
GLIMMIX macros in SAS System, but I would like to fit this kind of models 
in R. 

Could  anyone  help me on what package I should use to?

Thanks in advance

Francisco

[[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] Generalized Linear Mixed Models

2005-06-29 Thread Berton Gunter
submit the following in R:

RSiteSearch('glmm',restr='functions')

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
The business of the statistician is to catalyze the scientific learning
process.  - George E. P. Box
 
 

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of 
 [EMAIL PROTECTED]
 Sent: Wednesday, June 29, 2005 12:05 PM
 To: r-help@stat.math.ethz.ch
 Subject: [R] Generalized Linear Mixed Models
 
 Hello! 
 
 I am trying to fit a Generalized Linear Mixed Model, ordinally I use 
 GLIMMIX macros in SAS System, but I would like to fit this 
 kind of models 
 in R. 
 
 Could  anyone  help me on what package I should use to?
 
 Thanks in advance
 
 Francisco
 
   [[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


Re: [R] Generalized Linear Mixed Models

2005-06-29 Thread Ruben Roa
 -Original Message-
 From: [EMAIL PROTECTED]
 [mailto:[EMAIL PROTECTED] Behalf Of
 [EMAIL PROTECTED]
 Sent: 29 June 2005 18:05
 To: r-help@stat.math.ethz.ch
 Subject: [R] Generalized Linear Mixed Models
 
 
 Hello! 
 
 I am trying to fit a Generalized Linear Mixed Model, ordinally I use 
 GLIMMIX macros in SAS System, but I would like to fit this 
 kind of models in R. 
 
 Could  anyone  help me on what package I should use to?
 
 Thanks in advance
 
 Francisco

Si el modelo es de cualquier clase, revisa el package lme4.
Si el modelo es espacial en las familias Poisson o binomial, revisa package 
geoRglm.
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


[R] return NA

2005-06-29 Thread dscully




A-c(1,2,NA,7,5)
B-c(3,4,1,4,1)
C-c(6,5,6,NA,9)
D-c(8,7,4,6,2)
df1-cbind(A,B,C,D)
for(i in seq(1,ncol(df1)-1, by=2)) {
   ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] )  }


Tried several variations but none worked.  I wish to find any NA's in
column's 1 or 3 and change the numerical value to the right  of  the 
NA's .  In this case I wish to replace the 1 and 6 respectively with NA.
Any help would be appreciated.

__
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] return NA

2005-06-29 Thread Liaw, Andy
You need to use is.na(df1[,i]) to test for NA.

Andy

 From: [EMAIL PROTECTED]
 
 
 
 
 
 A-c(1,2,NA,7,5)
 B-c(3,4,1,4,1)
 C-c(6,5,6,NA,9)
 D-c(8,7,4,6,2)
 df1-cbind(A,B,C,D)
 for(i in seq(1,ncol(df1)-1, by=2)) {
ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] )  }
 
 
 Tried several variations but none worked.  I wish to find any NA's in
 column's 1 or 3 and change the numerical value to the right  of  the 
 NA's .  In this case I wish to replace the 1 and 6 
 respectively with NA.
 Any help would be appreciated.
 
 __
 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] return NA

2005-06-29 Thread Sundar Dorai-Raj


[EMAIL PROTECTED] wrote:
 
 
 
 A-c(1,2,NA,7,5)
 B-c(3,4,1,4,1)
 C-c(6,5,6,NA,9)
 D-c(8,7,4,6,2)
 df1-cbind(A,B,C,D)
 for(i in seq(1,ncol(df1)-1, by=2)) {
ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] )  }
 
 
 Tried several variations but none worked.  I wish to find any NA's in
 column's 1 or 3 and change the numerical value to the right  of  the 
 NA's .  In this case I wish to replace the 1 and 6 respectively with NA.
 Any help would be appreciated.
 

I think you want ?is.na.

for(i in seq(1, ncol(df1) - 1, by = 2))
   df1[, i + 1] - ifelse(is.na(df1[, i]), df1[, i], df1[, i + 1])

--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] predicted survival curve from a cox model

2005-06-29 Thread Lisa Wang
Hi there,

I have a predictor varible class which is a categorical variable and  
a ' coxph' is used to find the coeffients. How can I plot the predicted 
survival proportion based on this model?

Thanks

Lisa Wang
Princess Margaret Hospital
Toronto
tel 416 946 4501

__
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] plot (log scale on y-axis)

2005-06-29 Thread Jing Shen
I am planning to plot my data on log scale (y-axis). There is a
parameter in plot function, which is
plot( ..., log=y, ...)
While, the problem is that it is with base of e. Is there a way to let
me change it to 10 instead of e?

Thanks

__
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] return NA

2005-06-29 Thread Jim Brennan
Here is a way to do it without a loop that could save some time for a big
dataset.
 df1
  A B  C D
[1,]  1 3  6 8
[2,]  2 4  5 7
[3,] NA 1  6 4
[4,]  7 4 NA 6
[5,]  5 1  9 2

 df2-cbind(0,ifelse(is.na(df1),NA,0))[,-ncol(df1)-1]
 df2
A B  C
[1,] 0  0 0  0
[2,] 0  0 0  0
[3,] 0 NA 0  0
[4,] 0  0 0 NA
[5,] 0  0 0  0
 df3-df1+df2
 df3
  A  B  C  D
[1,]  1  3  6  8
[2,]  2  4  5  7
[3,] NA NA  6  4
[4,]  7  4 NA NA
[5,]  5  1  9  2

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of
[EMAIL PROTECTED]
Sent: June 29, 2005 3:50 PM
To: r-help@stat.math.ethz.ch
Subject: [R] return NA





A-c(1,2,NA,7,5)
B-c(3,4,1,4,1)
C-c(6,5,6,NA,9)
D-c(8,7,4,6,2)
df1-cbind(A,B,C,D)
for(i in seq(1,ncol(df1)-1, by=2)) {
   ifelse(df1[,i]==NA,df1[,i+1]==NA,df1[,] )  }


Tried several variations but none worked.  I wish to find any NA's in
column's 1 or 3 and change the numerical value to the right  of  the 
NA's .  In this case I wish to replace the 1 and 6 respectively with NA.
Any help would be appreciated.

__
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] sbrier (Brier score) and coxph

2005-06-29 Thread Stephen Henderson
Hello
I've decided to try and distill an earlier rather ill focused question to
try and elicit a response. Any help is greatly appreciated. Why does mod.cox
not work with sbrier whilst mod.km does? Can I make it work?



 data(DLBCL)
 DLBCL.surv-Surv(DLBCL$time,DLBCL$cens)
 
 mod.km-survfit(DLBCL.surv)
 mod.cox-survfit(coxph(DLBCL.surv~IPI, data=DLBCL))
 
 
 sbrier(DLBCL.surv, mod.km)
integrated Brier score 
 0.2076454 
attr(,time)
[1]   1.3 129.9
 sbrier(DLBCL.surv, mod.cox)
Error in switch(ptype, survfit = { : switch: EXPR must return a length 1
vector

Thanks in advance
Stephen Henderson



**
This email and any files transmitted with it are confidentia...{{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


Re: [R] sbrier (Brier score) and coxph

2005-06-29 Thread Prof Brian Ripley
Is this sbrier from package ipred?

The short answer is that it contains

 ptype - class(pred)

and assumes that is of length one.  For a survfit.coxph fit it is of class
c(survfit.cox, survfit).  I suspect from the help page that this is 
not supported, but you need to contact the authors (as the posting guide 
suggests).

On Wed, 29 Jun 2005, Stephen Henderson wrote:

 Hello
 I've decided to try and distill an earlier rather ill focused question to
 try and elicit a response. Any help is greatly appreciated. Why does mod.cox
 not work with sbrier whilst mod.km does? Can I make it work?



 data(DLBCL)
 DLBCL.surv-Surv(DLBCL$time,DLBCL$cens)

 mod.km-survfit(DLBCL.surv)
 mod.cox-survfit(coxph(DLBCL.surv~IPI, data=DLBCL))


 sbrier(DLBCL.surv, mod.km)
 integrated Brier score
 0.2076454
 attr(,time)
 [1]   1.3 129.9
 sbrier(DLBCL.surv, mod.cox)
 Error in switch(ptype, survfit = { : switch: EXPR must return a length 1
 vector

 Thanks in advance
 Stephen Henderson



 **
 This email and any files transmitted with it are confidentia...{{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


-- 
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 convert c:\a\b to c:/a/b

2005-06-29 Thread Spencer Graves
Thank You, Prof. Ripley!

  Both test1.R and test2.R worked for me just now, as did the 
following minor modification:

(x - readLines(stdin(), n=1))
D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt

  Thanks again.

  spencer graves

Prof Brian Ripley wrote:
 On Wed, 29 Jun 2005, David Duffy wrote:
 
 
I couldn't resist adding a more literal answer
 
 
 This can only work for escapes which are preserved.  The parser maps
 \n to a character (LF) and the deparser maps it back to \n.
 This happens to be true of \a \b \f \n \r \t \v \\ but no others.
 
 For example, \s is mapped to s, and there is no difference between \s and
 s in the parsed input.
 
 
unback - function(x) {
 chars - unlist(strsplit(deparse(x),))
 chars - chars[-c(1,length(chars))]
 paste(gsub(,/,chars),collapse=)
}

unback(\n)
 
 
unback(\s)
 
 [1] s
 
 Spencer Graves keeps on insisting there is a better way, but all the
 solutions are to avoid sending the string to the parser, and hence
 avoiding having the string directly in an R script.  This is common in 
 shell scripts, which use 'here' documents to avoid 'quoting hell'.
 
 We can do that in R too. Here are two variants I have not seen in the 
 thread
 
 test1.R:
 scan(, , allowEscapes=FALSE, n=1, quiet=TRUE)
 D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 catIx, \n, sep=)
 
 R --slave --vanilla  test1.R
 D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 
 (This one does not allow quoted strings.)
 
 test2.R:
 x - readLines(stdin(), n=1)
 D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 x - gsub('^(.*)$', \\1, x)
 cat(x, \n)
 
 R --slave --vanilla  test2.R
 D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 
 (This one allows surrounding quotes or not.)
 

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


Re: [R] How to convert c:\a\b to c:/a/b?

2005-06-29 Thread Gabor Grothendieck
One other comment.  Ninotech Path Copy, which can be found at:

  http://home.worldonline.dk/ninotech/

is a free Windows utility that appears in the Windows Explorer context
menu (i.e. it appears as the Copy Path menu entry when you right click 
any file in Windows Explorer).  I had forgotten about it since even though 
I had installed it a long time ago I hardly ever use it.  When one right clicks
and chooses Copy Path one is given a number of different ways in which
to copy the path.  Using the Setup submenu (i.e. right click any file and
choose Setup) one can define additional ways and its quite easy to define 
a new method which replaces backslashes with forward slashes as it 
copies -- in fact, there is a check box specifically for this.  

On 6/27/05, Spencer Graves [EMAIL PROTECTED] wrote:
 Hi, Gabor, James, et al.:
 
  Thanks for the help.  I'm unable to get scan('clipboard', what='',
 allowEscapes=FALSE)) to work reliably on my system using Rgui under
 Windows XP, R 2.1.1 patched.  However, file.choose() works like a charm.
 
  Thanks again, everyone.
  spencer graves
 
 Gabor Grothendieck wrote:
  Try using file.choose() to locate the file instead of Windows Explorer.  
  That
  will return the name in a form useable within R.
 
  On 6/27/05, Spencer Graves [EMAIL PROTECTED] wrote:
 
  Thanks, Dirk, Gabor, Eric:
 
  You all provided appropriate solutions for the stated problem.
 Sadly, I oversimplified the problem I was trying to solve:  I copy a
 character string giving a DOS path from MS Windows Explorer into an R
 script file, and I get something like the following:
 
  D:\spencerg\statmtds\R\Rnews
 
  I want to be able to use this in R with its non-R meaning, e.g., in
 readLine, count.fields, read.table, etc., after appending a file name.
 Your three solutions all work for my oversimplified toy example but are
 inadequate for the problem I really want to solve.
 
  Thanks,
  spencer graves
 
 Gabor Grothendieck wrote:
 
 
 On 6/27/05, Dirk Eddelbuettel [EMAIL PROTECTED] wrote:
 
 
 On 26 June 2005 at 20:30, Spencer Graves wrote:
 | How can one convert back slashes to forward slashes, e.g, 
 changing
 | c:\a\b to c:/a/b?  I tried the following:
 |
 |   gsub(, /, c:\a\b)
 | [1] c:\a\b
 
 This does work, provided you remember that single backslashed don't 
 exist
 as e.g. \a is a character in itself. So use doubles are you should be 
 fine:
 
 
 
 gsub(, /, c:\\a\\b)
 
 [1] c:/a/b
 
 
 
 Also, if one finds four backslashes confusing one can avoid the use
 of four via any of these:
 
 gsub([\\], /, c:\\a\\b)
 gsub(\\, /, c:\\a\\b, fixed = TRUE)
 chartr(\\, /, c:\\a\\b)
 
 __
 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
 
 
 
 
 --
 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


Re: [R] How to convert c:\a\b to c:/a/b

2005-06-29 Thread Gabor Grothendieck
Note that if you want to source it rather than than run it as a batch
job from the command line you will something like this.   The way
it works is that one puts the file name into comments that are
marked with tags and the script rereads itself as data picking out
the tagged lines and removing the tags from those lines.
script.name() returns the name of the script its running in and
my.stdin sets up a connection to the script as data after stripping out
all lines that do not match the tag and removing the tag from
those lines that do.  (This is a variation of some code that I previously
posted.)

# source the following from R

script.name - function() 
 showConnections()[as.character(eval.parent(quote(file), n=3)), description]
my.stdin - function( tag = ^#, script.name. = script.name())
 textConnection( sub(tag, , grep(tag,readLines(script.name.),value=TRUE)) )

x - readLines( my.stdin() )
#c:\a\b
cat(x, \n)



On 6/29/05, Spencer Graves [EMAIL PROTECTED] wrote:
 Thank You, Prof. Ripley!
 
  Both test1.R and test2.R worked for me just now, as did the
 following minor modification:
 
 (x - readLines(stdin(), n=1))
 D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 
  Thanks again.
 
  spencer graves
 
 Prof Brian Ripley wrote:
  On Wed, 29 Jun 2005, David Duffy wrote:
 
 
 I couldn't resist adding a more literal answer
 
 
  This can only work for escapes which are preserved.  The parser maps
  \n to a character (LF) and the deparser maps it back to \n.
  This happens to be true of \a \b \f \n \r \t \v \\ but no others.
 
  For example, \s is mapped to s, and there is no difference between \s and
  s in the parsed input.
 
 
 unback - function(x) {
  chars - unlist(strsplit(deparse(x),))
  chars - chars[-c(1,length(chars))]
  paste(gsub(,/,chars),collapse=)
 }
 
 unback(\n)
 
 
 unback(\s)
 
  [1] s
 
  Spencer Graves keeps on insisting there is a better way, but all the
  solutions are to avoid sending the string to the parser, and hence
  avoiding having the string directly in an R script.  This is common in
  shell scripts, which use 'here' documents to avoid 'quoting hell'.
 
  We can do that in R too. Here are two variants I have not seen in the
  thread
 
  test1.R:
  scan(, , allowEscapes=FALSE, n=1, quiet=TRUE)
  D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
  catIx, \n, sep=)
 
  R --slave --vanilla  test1.R
  D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 
  (This one does not allow quoted strings.)
 
  test2.R:
  x - readLines(stdin(), n=1)
  D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
  x - gsub('^(.*)$', \\1, x)
  cat(x, \n)
 
  R --slave --vanilla  test2.R
  D:\spencerg\dataPOWER\stats\Tukey\Boxplot_missing_Tukey2.txt
 
  (This one allows surrounding quotes or not.)
 
 
 --
 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-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] Extract fixed effects SE from lmer

2005-06-29 Thread La Sorte, Frank A.
Hi,
 
Does anyone know how to extract fixed effects SE values from generalized linear 
mixed models estimated using the lmer function in the lme4 library?  I searched 
attributes and structure with no luck.
 
Thanks
 
Frank A. La Sorte, Ph.D.
Department of Fisheries and Wildlife Sciences
University of Missouri
Columbia, MO 65211 USA

__
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] Fwd: Extract fixed effects SE from lmer

2005-06-29 Thread Douglas Bates
I forgot to cc: the list on this reply.

-- Forwarded message --
From: Douglas Bates [EMAIL PROTECTED]
Date: Jun 29, 2005 6:28 PM
Subject: Re: [R] Extract fixed effects SE from lmer
To: La Sorte, Frank A. [EMAIL PROTECTED]


On 6/29/05, La Sorte, Frank A. [EMAIL PROTECTED] wrote:
 Hi,

 Does anyone know how to extract fixed effects SE values from generalized 
 linear mixed models estimated using the lmer function in the lme4 library?  I 
 searched attributes and structure with no luck.


Try vcov(myLmerModel).  That returns a variance-covariance matrix for
the fixed effects.  The usual manipulations are used to get the
standard errors.  For the binomial and Poisson families you should use
the optional argument useScale=FALSE as these families do not have a
separately estimated scale parameter in the variance function.

In the show method for the summary.lmer class the relevant piece of code is

  corF - as(as(vcov(object, useScale = useScale), pdmatrix),
 corrmatrix)

__
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] deal package

2005-06-29 Thread Weiwei Shi
Hi,
I am wondering if anyone here used deal package in R to do the
bayesian network. I am curious about its scalability: how many
variables and how many observations can it handle in a reasonable
time. If you have some good experience, please share your data
configurations.


thanks,


-- 
Weiwei Shi, Ph.D

Did you always know?
No, I did not. But I believed...
---Matrix III

__
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] Finding out collinearity in regression

2005-06-29 Thread Young Cho
Hi, I am trying to find out a collinearity in
explanatory variables with alias(). 

I creat a dataframe:

dat - ds[,sapply(ds,nlevels)=2]
dat$Y - Response

Explanatory variables are factor and response is
continuous random variable. When I run a regression, I
have the following error:

fit - aov( Y ~ . , data = dat)
Error in contrasts-(`*tmp*`, value =
contr.treatment) :
contrasts can be applied only to factors with
2 or more levels

I think there is a dependency in explanatory
variables. So, I wanted to use alias to find out a
dependency in design matrix but I can't because I
cannot create fit in the first place.

One of examples I found is:

carprice1.lm - lm(gpm100 ~
Type+Min.Price+Price+Max.Price+Range.Price,data=carprice)
alias(carprice1.lm)

But, what if I can create lm object ? Then is there a
way to find out a dependency in design matrix? Thanks
a lot for help in advance!

-Young.

__
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] Finding out collinearity in regression

2005-06-29 Thread Kevin Wang
Hi,

Young Cho wrote:

 fit - aov( Y ~ . , data = dat)
 Error in contrasts-(`*tmp*`, value =
 contr.treatment) :
 contrasts can be applied only to factors with
 2 or more levels
 
 I think there is a dependency in explanatory
 variables. So, I wanted to use alias to find out a
 dependency in design matrix but I can't because I
 cannot create fit in the first place.

The error message actually looks like you have got (at least) a variable 
that only has 1 level, e.g. a factor with only one level.

Cheers,

Kev

-- 
Ko-Kang Kevin Wang
PhD Student
Centre for Bioinformation Science
Building 27, Room 1004
Mathematical Sciences Institute (MSI)
Australian National University
Canberra, ACT 2601
Australia

Homepage: http://wwwmaths.anu.edu.au/~wangk/
Ph (W): +61-2-6125-2431
Ph (H): +61-2-6125-7407
Ph (M): +61-40-451-8301

__
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] Finding out collinearity in regression

2005-06-29 Thread Simon Blomberg
At 01:45 PM 30/06/2005, Young Cho wrote:
Hi, I am trying to find out a collinearity in
explanatory variables with alias().

I creat a dataframe:

dat - ds[,sapply(ds,nlevels)=2]
dat$Y - Response

Explanatory variables are factor and response is
continuous random variable. When I run a regression, I
have the following error:

fit - aov( Y ~ . , data = dat)
Error in contrasts-(`*tmp*`, value =
contr.treatment) :
 contrasts can be applied only to factors with
2 or more levels

1. Sounds like at least one of your factors has only one level. This should 
be easy to spot.

2. Have you considered package perturb?

HTH,

Simon.



I think there is a dependency in explanatory
variables. So, I wanted to use alias to find out a
dependency in design matrix but I can't because I
cannot create fit in the first place.

One of examples I found is:

carprice1.lm - lm(gpm100 ~
Type+Min.Price+Price+Max.Price+Range.Price,data=carprice)
alias(carprice1.lm)

But, what if I can create lm object ? Then is there a
way to find out a dependency in design matrix? Thanks
a lot for help in advance!

-Young.

__
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

Simon Blomberg, B.Sc.(Hons.), Ph.D, M.App.Stat.
Centre for Resource and Environmental Studies
The Australian National University
Canberra ACT 0200
Australia
T: +61 2 6125 7800 email: Simon.Blomberg_at_anu.edu.au
F: +61 2 6125 0757
CRICOS Provider # 00120C

__
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] how to call egarch of sas in R

2005-06-29 Thread Nongluck Klibbua
I use R to generate data and I need to estimate the data by egarch (that
doesn't have in R). So how I can call egarch from SAS in R. 
Regards,
luck

__
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] moving correlation coef ?

2005-06-29 Thread vincent
Sundar Dorai-Raj a écrit :

 Perhaps ?cov.wt will work for you? Your example would be identical to:
 
 set.seed(1)
 X - rnorm(100); Y - rnorm(100)
 # using cov.wt
 rho1 - cov.wt(cbind(X, Y), 1:100, cor = TRUE)$cor[1, 2]
 # your weighting scheme
 rho2 - cor(X[rep(1:100, 1:100)], Y[rep(1:100, 1:100)])
 all.equal(rho1, rho2)
 # [1] TRUE
 
 HTH,
 
 --sundar

Thank you very much Sundar for your advices.
I will try this way.
Have a good day.
Vincent

__
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] Dispersion parameter in Neg Bin GLM

2005-06-29 Thread Edward McNeil
Hi,
Can someone tell me if it is possible to set the dispersion parameter constant 
when fitting a negative binomial glm in R? I've looked at the documentation and 
can't find the appropriate argument to pass.
In STATA I can type: nbreg depvar [indepvar...], offset(offset) 
dispersion(constant).
Thank you

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