Re: [R] Computing time for matrix addition or subtraction

2006-11-13 Thread Prof Brian Ripley
On Sun, 12 Nov 2006, YONGWAN CHUN wrote:

 I wonder by chance if there is a way to reduce computing time for matrix 
 addition or subtraction. With a lot of iterations, it would be helpful 
 to reduce a little amount time.

Yes, by making use of an optimized BLAS: see the R-admin manual.  On my 
(dual CPU) system this reduced your example from 36 to 6 seconds.

BTW, it is the calculation of PP that is taking the most of time, not as 
in your subject line.

 Simple example is as below

 n - 2000
 P - matrix(rnorm(n*n),n,n)
 PP - P %*% P
 M - diag(n) - P
 R - M + t(M) - diag(n) +  PP

 I would like to reduce time in calculating R.

-- 
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Computing time for matrix addition or subtraction

2006-11-13 Thread Matthias Kohl
not concerning your subject line, but function crossprod may be useful, too

Matthias

- original message 

Subject: Re: [R] Computing time for matrix addition or subtraction
Sent: Mon, 13 Nov 2006
From: Prof Brian Ripley[EMAIL PROTECTED]

 On Sun, 12 Nov 2006, YONGWAN CHUN wrote:
 
  I wonder by chance if there is a way to reduce computing time for matrix 
  addition or subtraction. With a lot of iterations, it would be helpful 
  to reduce a little amount time.
 
 Yes, by making use of an optimized BLAS: see the R-admin manual.  On my 
 (dual CPU) system this reduced your example from 36 to 6 seconds.
 
 BTW, it is the calculation of PP that is taking the most of time, not as 
 in your subject line.
 
  Simple example is as below
 
  n - 2000
  P - matrix(rnorm(n*n),n,n)
  PP - P %*% P
  M - diag(n) - P
  R - M + t(M) - diag(n) +  PP
 
  I would like to reduce time in calculating R.
 
 -- 
 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
 and provide commented, minimal, self-contained, reproducible code.
 

--- original message end 

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Need help in waveslim package: imodwt and universal.thresh.modwt

2006-11-13 Thread Airon Yiu
Hi:
  I have encountered problems with imodwt and universal.thresh.modwt and cannot 
find any reference in R Search.   Hope someone can give me some ideas:
   
  Starting with 
  modwt.la8 - modwt(xdata, la8, n.level=6)  -- this seems to work fine 
   
  (1)  ydata - imodwt(modwt.la8)
  will always give ydata as numeric(0) (no values) instead of being a time 
series data with the same lenght as xdata.  
   
  (2)   thred.la8 - universal.thresh.modwt(modwt.la8, max.level=4, hard=F) 
will give me error at:
  abs(wc.fine)  - cannot contain non-numeric field.
   
  Any help is much appreciated.
   
  Regards

 ___
 YM - Â÷½u°T®§
 
´Nºâ§A¨S¦³¤Wºô¡A§AªºªB¤Í¤´¥i¥H¯d¤U°T®§µ¹§A¡A·í§A¤Wºô®É´N¯à¥ß§Y¬Ý¨ì¡A¥ô¦ó»¡¸Ü³£ÉN¨«¥¢¡C

[[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
and provide commented, minimal, self-contained, reproducible code.


[R] doubt on Tukey HSD

2006-11-13 Thread Jose Kaitharan
dear all

effect of a, b and c on d, total 48 comparisons, got
one anova result in model1=aov(d~A*B*C). can we get
all the result in one command? or can we interpret the
whole comparisons from this result? how it work in
Tukey HSD?

jose




 


Access over 1 million songs.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Problems with metaMDS from vegan

2006-11-13 Thread Peter Roosen
Hello Spencer,

   1.  Have you tried to develop an example so small and simple that 
 it can be provided with a few lines of code in an email like this (as 
 suggested in the trailer to every email distributed by r-help)?

I did, but due to the simplicity of the example in my former posting you 
may have overlooked it.

 You should not criticize the maintainer 
 for not replying to your earlier email:  I did not see in your email any 
 information that would allow anyone to solve your problem unless s/he 
 just happened to be unbelievably lucky.

Well, the mentioning of my post to the maintainer was not intended as 
criticism, but as justification for the public posting itself.

 Please include 'sessionInfo()'

This does not reveal any additional information compared to that given 
in my former post, but anyhow (I noticed a new version of Vegan on 
Jari's HP and updated accordingly in the hope that it might cure the 
effect - which it didn't):


On the machine where the IDENTICAL CODE AND DATA is defunct:

  sessionInfo()
R version 2.2.1, 2005-12-20, i486-pc-linux-gnu

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

other attached packages:
 MASSvegan
7.2-24  1.9-3


On the machine where the code and data runs smoothly:

  sessionInfo()
R version 2.2.1, 2005-12-20, i486-pc-linux-gnu

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

other attached packages:
 MASSvegan
7.2-24 1.6-10


(The third machine [Win XP with R 2.4.0 and Vegan 1.8.3-x], show the 
defunct symptoms, is not available for testing right now)


   2.  Have you tried 'traceback()' after your failed call to 
 'metaMDS'?

traceback() on the defunct installation yields:

  traceback()
6: stop('x' must be an array of at least two dimensions)
5: rowSums(x, na.rm = TRUE)
4: any(rowSums(x, na.rm = TRUE) == 0)
3: vegdist(comm, method = distance, ...)
2: metaMDSdist(comm, distance = distance, autotransform = autotransform,
noshare = noshare, trace = trace, commname = commname, ...)
1: metaMDS(distab.dist, distance = bray, k, trymax = 50)


 With luck, this will identify the function actually 
 generating the error.  Then you can list that function and try find 
 rowSums(x, a.rm = TRUE).  Then hunt for x in that code.

That rowSums() code is identical on both installations, but the higher 
layers of invocation differ in differing amounts. I checked with mgdiff.

 The 
 function that includes this call to 'rowSums' expects 'x' to be an array 
 of at least two dimensions.  Either you pass it something that does NOT 
 meet that description or someplace x might be created from selecting 
 rows or columns of a larger array and either 0 or 1 rows were selected, 
 thereby reducing the dimension of the array to 0 or 1.

Yes, I guessed such a dependency already, but do not see the exact place 
where to dig into. The major code differences are in the vegdist() and 
the metaMDSdist() functions but I cannot point my finger to a certain 
code part where I estimate it to break - sorry!

   3.  If 1 and 2 fail to provide enlightenment, I would then try 
 the 'debug' functions (provided it was important enough for me to do so, 
 of course). ...

This is my next step to try as the first two do not seem to yield 
further insights for me, but I just wanted to return a timely reply to 
the other points you mentioned in your post and with a certain hope that 
it might give some longer term R users some ideas to communicate some ideas.

Thanks for your involvement,

Peter

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] A printing macro

2006-11-13 Thread Murray Jorgensen
I am exploring the result of clustering a large multivariate data set 
into a number of groups, represented, say, by a factor G.

I wrote a function to see how categorical variables vary between groups:

   ddisp - function(dvar) {
+  csqt - chisq.test(G,dvar)
+  print(csqt$statistic)
+  print(csqt$observed)
+  print(round(csqt$expected))
+  round(csqt$residuals)
+  }
 
   x - ceiling(4*runif(100))
   G - gl(4,1,100)
   ddisp(x)
X-squared
  6.235645
dvar
G1  2  3  4
   1 10  5  5  5
   2  6  9  5  5
   3  8  6  5  6
   4  7  4  4 10
dvar
G   1 2 3 4
   1 8 6 5 6
   2 8 6 5 6
   3 8 6 5 6
   4 8 6 5 6
dvar
G1  2  3  4
   1  1  0  0 -1
   2 -1  1  0 -1
   3  0  0  0  0
   4  0 -1  0  1
Warning message:
Chi-squared approximation may be incorrect in: chisq.test(G, dvar)

As I need to apply this function to a large number of variables x it 
would be helpful if the function printed x rather than the formal 
argument dvar. I have a vague idea that things like deparse() and 
substitute() will come into the solution but I have not yet come up with 
the right incantation. Any help appreciated!

Murray Jorgensen

-- 
Dr Murray Jorgensen  http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: [EMAIL PROTECTED]Fax 7 838 4155
Phone  +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] stepAIC for overdispersed Poisson

2006-11-13 Thread Murray Jorgensen
I am wondering if stepAIC in the MASS library may be used for model 
selection in an overdispersed Poisson situation. What I thought of doing 
was to get an estimate of the overdispersion parameter phi from fitting 
a model with all or most of the available predictors (we have a large 
number of observations so this should not be problematical) and then use 
stepAIC with scale = phi. Should this be OK?

Murray Jorgensen
-- 
Dr Murray Jorgensen  http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: [EMAIL PROTECTED]Fax 7 838 4155
Phone  +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Problems with metaMDS: Resolution!

2006-11-13 Thread Peter Roosen
The problem is solved:

With the help of the Vegan library maintainer the error was found 
immediately. It was a misinterpretation/overlooking of the provided 
documentation.

The problem lies here:

   meta - metaMDS(distab.dist, distance=bray, k, trymax=50)

The provided primary data matrix - distab.dist - was in my case a 
distance matrix, not a original data set matrix as demanded in the 
documentation.

As Jari Oksanen told me, several other users seem to have been trapped 
by the same pragmatic error also, which are reported in his newer 
versions of the Vegan library only. At least version 1.6-10 ran smoothly 
without reporting an error, thus leaving the user with the faulty notion 
of a correctly running program.

It seems that the pragmatic error of not only me stemmed from the fact 
that the frequently used methods cmdscale and isoMDS from the original R 
distribution scope require indeed *distance matrices* as their first 
input value. So, switching from them to the Vegan metaMDS, makes one 
easily overlook the important change in required data structure. The 
error condition that the newer versions of Vegan now returns helps to 
detect this faulty data provision.

Kind regards and many thanks to the list and the maintainer,

Peter

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] legend positioning problems

2006-11-13 Thread Schweitzer, Markus
Hello,

I have several Barplots I want to plot. I also want to include a legend.
Since my Barplots are very different I ve dicided to put the legend in
the top left corner.

Unfortunately, sometimes there is a part of a bar just below the legend.

This makes it difficult to see the legend itself or the barplot.

All I found out so far, is to make the backgroundcolor of the legend
white or so, but I would couldn t see the bar anymore.

Is there a way, to position the legend excel-like on the outside the
plot on the side? or
Is there a possibility to use sort of an alpha-factor to make the
background white but slightly invisible so that there is a contrast?

I hope you can help me out and thank you in advance!

#Example:
barplot(VADeaths)
legend(topleft, c(1x, 2x, 3x, 4x, 5x),
bg=white))




Markus Schweitzer - http://www.pokertips.tk

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] stepAIC for overdispersed Poisson

2006-11-13 Thread Prof Brian Ripley
On Mon, 13 Nov 2006, Murray Jorgensen wrote:

 I am wondering if stepAIC in the MASS library may be used for model
 selection in an overdispersed Poisson situation. What I thought of doing
 was to get an estimate of the overdispersion parameter phi from fitting
 a model with all or most of the available predictors (we have a large
 number of observations so this should not be problematical) and then use
 stepAIC with scale = phi. Should this be OK?

Well no, as that quasi-Poisson model does not have an AIC.

Remember AIC assumes maximum likelihood fitting, and you don't have a 
likelihood here (even for fixed phi).  The problem is that an 
'overdispersed Poisson situation' is a situation, not a model (in the 
usual sense of a probability measure over possible outcomes).

You could use a negative binomial GLM (with fixed theta).

-- 
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
and provide commented, minimal, self-contained, reproducible code.


[R] change \baselinestretch in Soutput ?

2006-11-13 Thread Christian Hoffmann
Hi there,

Is it possible to change \baselinestretch in Soutput ? What should

\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\footnotesize}

look like, if it is possible to effect the change here?

Google did not help here.

best
Christian
-- 
Dr. Christian W. Hoffmann,
Swiss Federal Research Institute WSL
Zuercherstrasse 111, CH-8903 Birmensdorf, Switzerland
Tel +41-44-7392-277 (office), -111(exchange), -215  (fax)
[EMAIL PROTECTED],  www.wsl.ch/staff/christian.hoffmann

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Need help in waveslim package: imodwt and universal.thresh.modwt

2006-11-13 Thread rdporto1
The first two commands worked fine to me. I used
xdata = rnorm(128).

To let us help you more, specify your R version,
some data, OS etc, as usually asked.

Rogerio Porto.

-- Cabeçalho original ---

De: [EMAIL PROTECTED]
Para: r-help@stat.math.ethz.ch
Cópia: 
Data: Sun, 12 Nov 2006 16:02:45 +0800 (CST)
Assunto: [R] Need help in waveslim package: imodwt and universal.thresh.modwt

 Hi:
   I have encountered problems with imodwt and universal.thresh.modwt and 
 cannot find any reference in R Search.   Hope someone can give me some ideas:

   Starting with 
   modwt.la8 - modwt(xdata, la8, n.level=6)  -- this seems to work fine 

   (1)  ydata - imodwt(modwt.la8)
   will always give ydata as numeric(0) (no values) instead of being a time 
 series data with the same lenght as xdata.  

   (2)   thred.la8 - universal.thresh.modwt(modwt.la8, max.level=4, hard=F) 
 will give me error at:
   abs(wc.fine)  - cannot contain non-numeric field.

   Any help is much appreciated.

   Regards
 
  ___
  YM - Â÷½u°T®§
  
 ´Nºâ§A¨S¦³¤Wºô¡A§AªºªB¤Í¤´¥i¥H¯d¤U°T®§µ¹§A¡A·í§A¤Wºô®É´N¯à¥ß§Y¬Ý¨ì¡A¥ô¦ó»¡¸Ü³£ÉN¨«¥¢¡C
 
   [[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
and provide commented, minimal, self-contained, reproducible code.


[R] Heatmap.2 in gplots

2006-11-13 Thread Jean Vidal
I used the function heatmap.2 some times ago, and it worked as intended (for
me...). I tried to reuse my old R program and I encounter some trouble. I
suppose that something was changed in newer versions.

One of the examples given in the help page shows the same problem :

-- this part wa
 ## For variable clustering, rather use distance based on cor():
data(USJudgeRatings)
symnum( cU - cor(USJudgeRatings) )

hU - heatmap.2(cU, Rowv=FALSE, symm=TRUE, col=topo.colors(16),
  distfun=function(c) as.dist(1 - c), trace=none)

## The Correlation matrix with same reordering:
hM - format(round(cU[hU[[1]], hU[[2]]], 2))
hM

# now with the correlation matrix on the plot itself

heatmap.2(cU, Rowv=FALSE, symm=TRUE, col=rev(heat.colors(16)),
 distfun=function(c) as.dist(1 - c), trace=none,
 cellnote=hM)

When running this correlation values are supposed to appear on the plot.
But, as one can see in this example, the same value is not associated with
the same color ( sometimes 1.00 is red or pale yellow).
Looks like some reordering done on colors was'nt done on cellnotes.

[[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
and provide commented, minimal, self-contained, reproducible code.


[R] GESCHAEFTSVORCHLAG

2006-11-13 Thread Ehret Olds
Sehr geehrter
 
Geschaefts-Vorschlag!

Sicher sind Sie verwundert, diese Nachricht von jemandem zu erhalten, den
sie nicht personlich kennen. Der Grund weshalb ich mich an Sie wende ist das
ich Ehret Olds bin, der erstgeborene Sohn von Martin Olds, einem der
bekanntesten schwarzen Farmer in Zimbabwe.
Er wurde vor kurzem im Laufe des Landstreites in meinem Land ermordet. Ich
erhielt Ihre Adresse durch das Internet und habe mich entschieden mich an
Sie zu wenden.

Bevor mein Vater starb, nahm er mich mit nach Johannesburg um dort einen
Betrag von 8,5 Millionen US$ (Acht Millionen funfhunderttausend US Dollar)
bei einer privaten Sicherheitsfirma zu deponieren. Er hat die drohende
Gefahr geahnt und daher diesen Betrag in einer Kassette dort hinterlegt, die
er als normales Wertgut deklarierte um keinen Verdacht zu erwecken und
ausserdem die Kosten fuer die Hinterlegung gering zu halten.

Niemand weis also von dem Bargeld in dieser Kassette. Ursprunglich war
dieser Betrag fuer neue Maschinen und Chemikalien für seine Farm gedacht,
ausserdem wollte er neues Farmland in Schwerz erwerben und erschliesen.
Das Landproblem entstand als der Präsident von Zimbabwe, Mr. Robert Mugabe,
seine neue Landreform durchsetzte bei der nur die reichen weisen Farmer und
einige wenige schwarze Farmer Vorteile verschafft wurden, fuer alle anderen
war die Existenz nun bedroht. Daraus resultierten dann schlielich die
Unruhen die sich ausbreiteten und denen viele Menschen zum Opfer fielen und
ihr Leben lassen mussten. Leider ist mein Vater ebenfalls eines dieser Opfer
das sein Leben verlor. Vor diesem Hintergrund bin ich mit der Familie unter
lebensgefährlichen Umständen aus Zimbabwe nach Holland geflohen, hier haben
wir politisches Asyl beantragt. Das in Südafrika zurückgelassene Geld wollen
wir nun der Sicherheit wegen auf ein Konto hier übertragen.
Das Gesetz in Holland verbietet es Menschen die politisches Asyl beantragt
haben während dieses Verfahrens ein Konto zu eröffnen oder
irgendwelche Transfers/Transaktionen die Über die hollandische Grenze
hinausgehen abzuwickeln.
Als der erstgeborene Sohn bin ich nun für meine ganze Familie veranwortlich
und habe die Rolle meines Vaters als Bewacher und Beschitzer Übernommen.

Nun stehe ich vor dem Problem dieses Geld ohne Wissen der afrikanischen
Regierung hierher zu transferieren, ansonsten wurde man uns den gesamten
Betrag enteignen, dieses Geld ist aber alles das uns noch geblieben ist. Die
Südafrikanische Regierung unterstützt wohl die Regierung Zimbabwes, so dass
auch dort nichts von meinem Vorhaben bekannt werden darf.
Als Geschäftsmann suche ich nun nach einem Partner dem ich voll und ganz
vertrauen kann, dem ich meine und auch die Zukunft meiner ganzen Familie
anvertrauen kann. Ich möchte Sie wissen lassen, dass ein Plan völlig ohne
Risiko ist, falls
Sie mir und meiner Familie helfen wollen. Alles worum ich Sie bitte, ist
eine Vereinbarung mit der hollandischen Niederlassung der Sicherheitsfirma
zu treffen, das hinterlegte Wertgut auszuhändigen. Ich meinerseits habe
bereits Anweisung gegeben die Kassette von Südafrika nach Holland zu senden.
Vorher allerdings müssen noch einige wichtige Formalitäten erledigt werden,
wie die Änderung des Begunstigten an diesem Wertgut. Weiterhin habe ich vor
diese Geld nach Erhalt möglichst gut und gewinnbringend zu investieren. Ich
möchte Ihnen zwei Vorschläge unterbreiten.Erstens biete ich Ihnen einen Teil
der Gesamtsumme an, wenn Sie bereit wären, Ihr eigenes Konto für diese
Transaktion Verfügung zu stellen.Oder aber Sie sind an einer Partnerschaft
mit mir interessiert um diese Summe profitabel in Ihrem Land zu
investieren.Egal welcher Vorschlag Ihnen mehr zusagt, zögern Sie nicht mir
Ihre Entscheidung mitzuteilen.

Fuer alle entstehenden Unkosten habe ich 5% des Gesamtbetrages ingeplant. 
Sollten Sie nicht an einer Partnerschaft mit Mir interessieren sein, biete
ich Ihnen 25% des Betrages an und werde die verbleibenden 70% fuer meine
Investitionen in Ihrem Land nutzen.
Sie können mich jederzeit unter dieser e-Mail Adresse erreichen: 
( [EMAIL PROTECTED] ) Ich bitte Sie jedoch um Ihre absolute Verschwiegenheit
in dieser
Angelegenheit Dritten gegenüber. Gott schutze Sie!
 
Hochachtungsvoll,
Ihr ergebener,
Ehret Olds.

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Fetching Intraday data from Bloomberg

2006-11-13 Thread Shubha Vishwanath Karanth
Hi Everyone.

 

I am downloading intraday Bloomberg data from R.

 

The code I give is:

 

library(zoo)

library(chron)

library(RBloomberg)

 

conn-blpConnect(show.days=trading,na.action=previous.days,periodici
ty=daily)

dat-blpGetData(conn, VG1 Index, c(LAST_PRICE),
start=as.chron(as.Date(2006-9-01, %Y-%m-%d)),barfields=OPEN,
barsize=10,retval=data.frame)

blpDisconnect(conn)

write.table(dat,Z:\h1.csv,sep=,,row.names=F,col.names=T)

 

 

Now, though I give show.days=trading, in the above code, my intraday
downloads also have Saturdays and Sunday data (non-active days). Why is
this so? Is it not updated in R? Note that this happens only for
intraday downloads and not for the usual historical downloads where time
is not considered.

 

Please answer me ASAP...

 

Thanks in advance.

 

 


[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Confidence intervals for relative risk

2006-11-13 Thread Terry Therneau

Wolfgang,
  It is common to handle relative risk problems using Poisson regression.
  In your example you have 8 events out of 508 tries, and 0/500 in the second
data set.

 tdata - data.frame(y=c(8,0), n=c(508,500), group=1:0)
 fit - glm(y ~ group + offset(log(n)), data=tdata, family=poisson)

  Because of the zero, the standard beta/se(beta) confidence intervals don't
work.  In fact, the actual MLE for the ratio is infinity, which the above
glm decides is exp(33.85) = 5* 10^14 --- which is close enough to infinity
for me.  What you want is the lower limit of the interval, which you can
find with a profile likelihood.  That is, look at the curve of deviance vs
beta, draw a horizontal line 3.84 units down from the top (chisq on 1 df),
and where it itersects the curve is your confidence limit.

  For the above, since it is 2 groups with 2 coefficients, the deviance of
the full model happens to be 0.  Your approximation gave a lower limit of
.97 which is about exp(0), so I'll guess a solution between -2 and 2.

 xx - seq(-2, 2, length=21)
 for (i in 1:21) {
fit - glm(y ~ offset(group*xx[i] + log(n)), poisson, tdata)
print(c(xx[i], fit$deviance))
}

[1]  -2.0  33.80736
[1]  -1.8  31.02994
[1]  -1.6  28.33138
[1]  -1.4  25.72327
[1]  -1.2  23.21769
[1]  -1.0  20.82691
[1]  -0.8  18.56281
[1]  -0.6  16.43629
[1]  -0.4  14.45668
[1]  -0.2  12.63108
[1]  0.0 10.96387
[1] 0.2 9.45639
[1] 0.40 8.106805
[1] 0.60 6.910274
[1] 0.80 5.859303
[1] 1.00 4.944278
[1] 1.20 4.154088
[1] 1.40 3.476757
[1] 1.6 2.90003
[1] 1.8 2.41186
[1] 2.00 2.000785

 So the exact lower limit is somewhere between exp(1.4) and exp(1.2), which
is around 3.6.  One can easily refine this by tucking the fit into an
iterative search, or just resetting xx to a new range.

  The offset(log(n)) part of the model, BTW, is a well known trick in poisson
models.  It has to do with the fact that the likelihood is in terms of the
number of events, but we want coefficients in terms of rates, and E(y) = rate*n.
Adding an offset of xx[i] * group fits a model with the group effect fixed at
xx, but allowing the intercept to vary.  
  For more conceptual depth, you could look up 'rate regression' or 
'standardized mortality ratio' in the Encyclopedia of Biostatistics --- it is
in the latter computation that this approach is quite common.  The idea of
adding a fraction is found in Anscombe(1949), Transformations of Poisson, 
binomial and negative-binomial data.  Biometrika, vol 35, p246-254, but for
each single estimate not for the ratio.

Terry Therneau
Mayo Clinic

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] segfault in RMySQL dbConnect error handling

2006-11-13 Thread Alex Brown
Hi there.  I see in a post from 2002 that you got the following  
problem with RMySQL:
  con - dbConnect(m)

Process R segmentation fault at Wed Aug 28 08:21:11 2002

I have the same problem today:

drv=dbDriver(MySQL)

dbConnect(drv) # or with pretty much any other failing options

Program received signal SIGSEGV, Segmentation fault.
0x2b28d60c9fa0 in strlen () from /lib/libc.so.6
(gdb) bt
#0  0x2b28d60c9fa0 in strlen () from /lib/libc.so.6
#1  0x2b28d609c7f4 in vfprintf () from /lib/libc.so.6
#2  0x2b28d60b6a69 in vsprintf () from /lib/libc.so.6
#3  0x2b28d60a1f76 in sprintf () from /lib/libc.so.6
#4  0x2b28d7f66b1b in ?? ()
#5  0x00fbda01 in ?? ()
#6  0xa285 in ?? ()
#7  0x011eda68 in ?? ()
#8  0x00d982a0 in ?? ()
#9  0x in ?? ()

(gdb) info regi
rax0x8d 141
rbx0x1  1
rcx0x3  3
rdx0x7f9ccb08   140737481853704
rsi0x2b28d7f67f01   47454421942017 (argument 2)  
nrecognised R/S type for group
rdi0x8d 141 (argument 1)
rbp0x7f9cc990   0x7f9cc990
rsp0x7f9cc2f8   0x7f9cc2f8

(gdb) c

*** caught segfault ***
address 0x8d, cause 'memory not mapped'

Traceback:
1: .Call(RS_MySQL_newConnection, drvId, con.params, groups,  
default.file, PACKAGE = .MySQLPkgName)
2: mysqlNewConnection(drv, ...)
3: .class1(object)
4: .class1(object)
5: is(object, Cl)
6: .valueClassTest(standardGeneric(dbConnect), DBIConnection,  
dbConnect)
7: dbConnect(drv)

Possible actions:
1: abort (with core dump)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection:
It appears that it is crashing while handling an error relating to  
the database.

If I give it parameters which allow it to connect, it does not segfault.



System information:

XEON x86_64 with VT

Ubuntu 6.06.1 LTS \n \lRN)

MySQL 5.0.22-Debian_0ubuntu6.06.2-log

R 2.4.0 R version 2.4.0 (2006-10-03) from deb http://cran.R- 
project.org/bin/linux/ubuntu dapper/

It's a pristine 2.4.0 R install with the following libraries:

DBI_0.1-11.tar.gz  RMySQL_0.5-10.tar.gz

installed in a custom library location /export/downloads/R/packages

(normal location seems to behave the same)



note: it also failed on 2.3.0 on this machine



--Alex Brown

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] random forest regression

2006-11-13 Thread Naiara Pinto
Dear all,

 I am doing a regression in ramdomForest, using the option sampsize reduce
the number of records used to produce the randomForest object.
The manual says For classification, if sampsize is a vector of the length
the number of strata, then sampling is stratified by strata, and the
elements of sampsize indicate the numbers to be drawn from the strata.  I
need my sampling to be done with factors, but I am doing a regression. Does
anyone know a way to do that?

Thanks,

Naiara.

[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R and Fortran 9x -- advice

2006-11-13 Thread Mike Prager
Tamas K Papp [EMAIL PROTECTED] wrote:

 I found some bottlenecks in my R code with Rprof.  First I wanted to
 rewrite them in C, but a colleague keeps suggesting that I learn
 Fortran, so maybe this is the time to do it...
 
 1) I hear bad things about Fortran.  Sure, F77 looks archaic, but
 F90/95 seems nicer.  Is it worth learning, especially when I plan to
 use it mainly from R?  Dusting off my C knowledge would take a bit of
 time too, and I never liked C that much.

I'll answer this from the perspective of someone who uses
Fortran 95 regularly.  It is a modern language, far more
reliable and flexible than Fortran 77 and quite well suited to
most scientific problems.  I do think it's worth learning,
particularly if C is not to your taste.

Two free compilers for Fortran 95 are available.  It seems that
g95 is complete, while gfortran is nearing completion.  There
are also several high-quality commercial compilers, some of
which are free under certain operating systems and/or conditions
and others of which (Lahey) are quite inexpensive if one is
willing to work from the command line or one's own editor.

I can't address questions of R interoperability  -- not
something I've done.

-- 
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Can I execute the content of a character vector?

2006-11-13 Thread Luiz Rodrigo Tozzi
Hi

I want to know if there is any possibility of executing the content of a
vector, for example:

example=c(Test,1,0,0,0,seq(14,42,by=2),0,0,1)

i want to know if there is anything like execute(example[6])

i really need this because this object example is created from a parameter
file with names, flags and this seq somethimes will be something like
c(1,99,3) or even c()

executing the content would be an easier step for me

thanks!

-- 

 Luiz Rodrigo Lins Tozzi


[EMAIL PROTECTED]
  (21)91318150


http://luizrodrigotozzi.multiply.com/
http://www.flogao.com.br/luizrodrigotozzi

[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] [R-SIG-Finance] looking for functions that can test/estimate CAPM, APT, Fama's factor model, etc.

2006-11-13 Thread Michael
Hi Brian,

thanks a lot for your pointers! I've taken a look at the book and the R
example website. That's super! Some of the examples there are very good.

Yet I am still looking for Fama 3 factor model and Ross' APT implementation.
The concept is not hard per se, however I am not sure how to classify some
companies as H, and some companies as L and others as Value companies, and
the others as Growth companies. There are a lot of implementation details
that I am not sure of. Thus I guess learning from other people's
implementation is a better approach for me as a green-hand.

(The Rmetrics are not very complete I guess.)

Thanks a lot Brian and thanks everybody...





On 11/13/06, Brian G. Peterson [EMAIL PROTECTED] wrote:

 On Sunday 12 November 2006 22:41, Michael wrote:
   I am also looking for interesting statistical experiments about
   testing and estimating CAPM, APT, Fama models, etc.
   using R using financial series data...
   please give me some pointers... I have been searching
   the R archives for the past a few hours and I vaguely got to know
   that there are programs do these interesting statistical things, but
   I just could not find where are they...

 Look at at the 'portfolio' and 'RMetrics' packages.  RMetrics is actually
 a group of packages, not just a single module.  These should give you
 CAPM and more.

 I can also recommend this excellent introductory book:

 http://www.amazon.com/Statistics-Finance-Introduction-David-Ruppert/dp/0387202706
 Statistics and Finance: An Introduction by David Ruppert

 with R examples available here:
 http://www.stat.tamu.edu/~ljin/Finance/stat689-R.htm

 Once you've worked on this for a while, if you need pointers on a specific
 problem, this list may be able to be of greater assistance.

 Regards,

- Brian

 ___
 [EMAIL PROTECTED] mailing list
 https://stat.ethz.ch/mailman/listinfo/r-sig-finance


[[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
and provide commented, minimal, self-contained, reproducible code.


[R] Problem with fitted.value function

2006-11-13 Thread AC
Who can help me in finding the right fitted-value-function?Based upon a 
linear model (2-dimensional, 1 result-variable and 1 input-variable).
I found the perfect regression model (not always simply linear) by the 
function lm(). Based on the given input and output tuples. For these exact 
input values I can get the exact fitted values by the function 
fitted.values.
The problem is, once I want to extrapolate with the automaticly found 
function (not always the same) and new input-variables I don't know how to 
calculate automaticly the fitted values. I don't know how to use the fitted 
values
function with a given function and given input-variables but yet unknown 
result-values.

Who can help me get the right fitted value function and tell me how to use 
it.

Thanks in advance for any help.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Can I execute the content of a character vector?

2006-11-13 Thread Roger Bivand
On Mon, 13 Nov 2006, Luiz Rodrigo Tozzi wrote:

 Hi
 
 I want to know if there is any possibility of executing the content of a
 vector, for example:
 
 example=c(Test,1,0,0,0,seq(14,42,by=2),0,0,1)
 
 i want to know if there is anything like execute(example[6])

You can say:

eval(parse(text=example[6]))

but it is difficult to avoid shooting people in the feet with it, unless 
that is what you want ...

eval(parse(text=example[1]))


 
 i really need this because this object example is created from a parameter
 file with names, flags and this seq somethimes will be something like
 c(1,99,3) or even c()
 
 executing the content would be an easier step for me
 
 thanks!
 
 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: [EMAIL PROTECTED]

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


Re: [R] Can I execute the content of a character vector?

2006-11-13 Thread Gabor Grothendieck
On 11/13/06, Roger Bivand [EMAIL PROTECTED] wrote:
 On Mon, 13 Nov 2006, Luiz Rodrigo Tozzi wrote:

  Hi
 
  I want to know if there is any possibility of executing the content of a
  vector, for example:
 
  example=c(Test,1,0,0,0,seq(14,42,by=2),0,0,1)
 
  i want to know if there is anything like execute(example[6])

 You can say:

 eval(parse(text=example[6]))

 but it is difficult to avoid shooting people in the feet with it, unless
 that is what you want ...

 eval(parse(text=example[1]))


We could enclose strings to be evaluated in backticks:

example - c(Test,1,0,0,0,`seq(14,42,by=2)`,0,0,1)

pat - ^`(.*)`$
f - function(x)  # strip backticks in 1st and last posn and evaluate
   if (regexpr(pat, x)  0)
  eval.parent(parse(text = sub(pat, \\1, x)))
   else x

example - c(Test,1,0,0,0,`seq(14,42,by=2)`,0,0,1)
lapply(example, f)



 
  i really need this because this object example is created from a parameter
  file with names, flags and this seq somethimes will be something like
  c(1,99,3) or even c()
 
  executing the content would be an easier step for me
 
  thanks!
 
 

 --
 Roger Bivand
 Economic Geography Section, Department of Economics, Norwegian School of
 Economics and Business Administration, Helleveien 30, N-5045 Bergen,
 Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
 e-mail: [EMAIL PROTECTED]

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


__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R and Fortran 9x -- advice

2006-11-13 Thread Christos Hatzis
Mike,

Can you recommend any good books on Fortran 90/95?  I had been an old user
of Fortran 77 but haven't followed the developments in the last 15 years or
so...

Thanks.

Christos Hatzis, Ph.D.
Nuvera Biosciences, Inc.
400 West Cummings Park
Suite 5350
Woburn, MA 01801
Tel: 781-938-3830
www.nuverabio.com
 


-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Mike Prager
Sent: Monday, November 13, 2006 1:00 PM
To: r-help@stat.math.ethz.ch
Subject: Re: [R] R and Fortran 9x -- advice

Tamas K Papp [EMAIL PROTECTED] wrote:

 I found some bottlenecks in my R code with Rprof.  First I wanted to 
 rewrite them in C, but a colleague keeps suggesting that I learn 
 Fortran, so maybe this is the time to do it...
 
 1) I hear bad things about Fortran.  Sure, F77 looks archaic, but
 F90/95 seems nicer.  Is it worth learning, especially when I plan to 
 use it mainly from R?  Dusting off my C knowledge would take a bit of 
 time too, and I never liked C that much.

I'll answer this from the perspective of someone who uses Fortran 95
regularly.  It is a modern language, far more reliable and flexible than
Fortran 77 and quite well suited to most scientific problems.  I do think
it's worth learning, particularly if C is not to your taste.

Two free compilers for Fortran 95 are available.  It seems that
g95 is complete, while gfortran is nearing completion.  There are also
several high-quality commercial compilers, some of which are free under
certain operating systems and/or conditions and others of which (Lahey) are
quite inexpensive if one is willing to work from the command line or one's
own editor.

I can't address questions of R interoperability  -- not something I've done.

--
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.

__
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
and provide commented, minimal, self-contained, reproducible code.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R and Fortran 9x -- advice

2006-11-13 Thread Ravi Varadhan
Metcalf and Reid - FORTRAN 90/95 Explained

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]

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 





-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Christos Hatzis
Sent: Monday, November 13, 2006 2:05 PM
To: 'Mike Prager'; r-help@stat.math.ethz.ch
Subject: Re: [R] R and Fortran 9x -- advice

Mike,

Can you recommend any good books on Fortran 90/95?  I had been an old user
of Fortran 77 but haven't followed the developments in the last 15 years or
so...

Thanks.

Christos Hatzis, Ph.D.
Nuvera Biosciences, Inc.
400 West Cummings Park
Suite 5350
Woburn, MA 01801
Tel: 781-938-3830
www.nuverabio.com
 


-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Mike Prager
Sent: Monday, November 13, 2006 1:00 PM
To: r-help@stat.math.ethz.ch
Subject: Re: [R] R and Fortran 9x -- advice

Tamas K Papp [EMAIL PROTECTED] wrote:

 I found some bottlenecks in my R code with Rprof.  First I wanted to 
 rewrite them in C, but a colleague keeps suggesting that I learn 
 Fortran, so maybe this is the time to do it...
 
 1) I hear bad things about Fortran.  Sure, F77 looks archaic, but
 F90/95 seems nicer.  Is it worth learning, especially when I plan to 
 use it mainly from R?  Dusting off my C knowledge would take a bit of 
 time too, and I never liked C that much.

I'll answer this from the perspective of someone who uses Fortran 95
regularly.  It is a modern language, far more reliable and flexible than
Fortran 77 and quite well suited to most scientific problems.  I do think
it's worth learning, particularly if C is not to your taste.

Two free compilers for Fortran 95 are available.  It seems that
g95 is complete, while gfortran is nearing completion.  There are also
several high-quality commercial compilers, some of which are free under
certain operating systems and/or conditions and others of which (Lahey) are
quite inexpensive if one is willing to work from the command line or one's
own editor.

I can't address questions of R interoperability  -- not something I've done.

--
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.

__
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
and provide commented, minimal, self-contained, reproducible code.

__
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
and provide commented, minimal, self-contained, reproducible code.

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] 95% CI of spectrum

2006-11-13 Thread gwang
Dear R-list Members,

Sorry if I have double posted this message.  I am not sure if my previous 
message got through.

Can someone explain me what exactly is the vertical bar of the cross in the 
upper right corner of the spectrum plot produced by spec.plot?  I know it is 
the 95% confidence of the spectrum, but do not understand what it is exactly.

Based on the section 9.5 of Bloomfield 2000, I understand how the pointwise 
95% CI of spectrum at a certain frequency is computed.  Is the vertical bar of 
the cross for the mean spectrum or total spectrum?  What is the exact 
interpretation of the vertical bar?  Many thanks in advance!

Sincerely,

Guiming Wang

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] help with syntax of nlme call.

2006-11-13 Thread Bill Shipley
I am getting an error message in a call to nlme and cannot understand what
is happening.  I explain the steps below in the hope that someone can
explain the error and how to correct it.
 
STEP 1: Data set: name: marouane.data.  This is a data frame whose first few
lines are as follows:
 
 marouane.data[1:13,]
   species plant leaf irradiance photosynthesis chlorophyll
1  Asclepias.incarnata 11  0  -2.091359 0.02619
2  Asclepias.incarnata 11 50   1.153241 0.02619
3  Asclepias.incarnata 11100   2.241963 0.02619
4  Asclepias.incarnata 11200   3.541258 0.02619
5  Asclepias.incarnata 11300   5.012547 0.02619
6  Asclepias.incarnata 11400   5.710689 0.02619
7  Asclepias.incarnata 11600   6.632571 0.02619
8  Asclepias.incarnata 11800   7.314284 0.02619
9  Asclepias.incarnata 11   1000   8.092402 0.02619
10 Asclepias.incarnata 11   1310   8.110368 0.02619
11 Asclepias.incarnata 12  0  -2.051934 0.02707
12 Asclepias.incarnata 12 50   1.213854 0.02707
13 Asclepias.incarnata 12100   2.271453 0.02707
   absorbance nitrogen   sla
1  0.8816 18.35913 77.53
2  0.8816 18.35913 77.53
3  0.8816 18.35913 77.53
4  0.8816 18.35913 77.53
5  0.8816 18.35913 77.53
6  0.8816 18.35913 77.53
7  0.8816 18.35913 77.53
8  0.8816 18.35913 77.53
9  0.8816 18.35913 77.53
10 0.8816 18.35913 77.53
11 0.8813 18.35913 77.53
12 0.8813 18.35913 77.53
13 0.8813 18.35913 77.53

The nesting structure is species (25), plants within species (4 each
species) and leaves within plants (2 each plant).  There is only 1 missing
value in the data set.
 
STEP 2: I constructed a self starting function (photo) that gives the
nonlinear function and provides initial estimates of the three parameters.
This self starting function works.  I use this to call the nlsList function,
which fits separate nonlinear functions to each species (I am only working
on a 2 level model so far: between and within species):
 
fit-nlsList(photosynthesis~photo(irradiance,Am,Q,LCP)|species,
+ data=marouane.data,na.action=na.omit)

This works and I get the three parameter estimates per species.
 
STEP 3: use nlsList to call nlme.  I use the nlsList object (fit) to fit a
variance components model:
 
fit.nlme-nlme(fit)
 
This works (at least it runs without an error message).  To see what syntax
is used, I extract the call:
 
 fit.nlme$call
nlme.formula(model = photosynthesis ~ photo(irradiance, Am, Q, 
LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~ 
1), random = list(species = c(2.41887429868478, -1.53351210977838, 
2.47282942435836, 0, 0, 0)), groups = ~species, start = list(
fixed = c(11.7384877502532, 0.109091641284836, 9.91905527125462
)), na.action = na.omit)

 
Question 1: the syntax to the random argument seems wrong!  It should be
something like: random=(list(Am~1,Q~1,LCP~1)).  I have no idea where the 6
numbers (2.41887...) are coming from in the random argument.  Can someone
explain how the nlme function obtains such an argument for random?  The
values in fixed are the estimated fixed effects from nlsList.
 
If I follow the instructions in Pinheiro  Bates, the call should be (with a
self-starting function):
nlme(model = photosynthesis ~ photo(irradiance, Am, Q, 
LCP), data = marouane.data, fixed = list(Am ~ 1, Q ~ 1, LCP ~ 
1), random = list(Am~1,Q~1,LCP~1),groups=~species,
na.action=na.omit).
When I use this, I get an error message:
 
  nlme(model=photosynthesis~photo(irradiance,Am,Q,LCP),data=marouane.data,
+  fixed=list(Am~1,Q~1,LCP~1),groups=~species,
+ random=list(Am~1,Q~1,LCP~1),
+  na.action = na.omit)
Error in nlmeCall[[i]] - NULL : subscript out of bounds

 
Bill Shipley
 

[[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
and provide commented, minimal, self-contained, reproducible code.


[R] indexing question

2006-11-13 Thread Leeds, Mark \(IED\)
I have the following set of indices, call it idx, that correspond to the
indices of a vector say temp.

 [1]   31   36   41   61   66   71   91   96  101  121  126  131  151
156  161  181  186  191  211  216  221  241  246  251  271  276  281
301  306  311  331  336  341  361  366
 [36]  371  391  396  401  421  426  431  451  456  461  481  486  491
511  516  521  541  546  551  571  576  581  601  606  611  631  636
641  661  666  671  691  696  701  721
 [71]  726  731  751  756  761  781  786  791  811  816  821  841  846
851  871  876  881  901  906  911  931  936  941  961  966  971  991
996 1001 1021 1026 1031 1051 1056 1061
[106] 1081 1086 1091  1116 1121 1141 1146 1151 1171 1176 1181 1201
1206 1211 1231 1236 1241 1261 1266 1271 1291 1296 1301 1321 1326 1331
1351 1356 1361 1381 1386 1391 1411 1416
[141] 1421


I want to calculate temp[36] - temp[31] and temp[41] - temp[36]

Similarly, temp[66] - temp[61] and temp[71] - temp[66]
.
.
.
.
Similarly temp[1416]-temp[1411]
  temp[1421] - temp[1416]


I'm doing this because the above subractions represent pairs of returns
and the correlations between them wil be calculated eventually.

In other words, eventually I will have

X_36_31  ( i.e : temp[36] - temp[31] )
X_66-61
X_96-91
.
.
.
.
.
.
.
X_1411-1416 

as one vector and

Y_41-36
Y_71-66
Y_101-96
.
.
.
.
.
Y_1416_1421

as another vector.

and will calculate the correlation between the two vectors in order to
get one number.


The point is I am  really only using the indices 31, 61, 91 etc as
anchor's so a regular diff(temp[idx]) won't work because it will diff
all 
the elements that are next to each other ? This is a weird problem. I'm
still thinking about it. I'm hoping to figure it out before someone
sends me something but I won't mind so much if I get an external
solution first. I have no pride.


This is not an offer (or solicitation of an offer) to buy/se...{{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
and provide commented, minimal, self-contained, reproducible code.


[R] hybrid in fisher.test broken?

2006-11-13 Thread Jonathan Dushoff
The hybrid feature in fisher.test looks to me like an excellent way to
analyze my two-way tables.   The only problem is that it does not seem
to be implemented.  Am I right about this?

An example is pasted below.  I note that I get the warning message only
when I shouldn't: for a 2x2 table hybrid seems to be ignored without
warning.  In no case does fisher.test seem to be checking Cochran
criteria as promised.

Any help will be appreciated.

JD

--


R : Copyright 2006, The R Foundation for Statistical Computing
Version 2.3.1 (2006-06-01)

.
.
.

 m = matrix(10*1:9, nc=3)
 fisher.test(m, hybrid=TRUE)
Error in fisher.test(m, hybrid = TRUE) : FEXACT error 6.
LDKEY is too small for this problem.
Try increasing the size of the workspace.
In addition: Warning message:
'hybrid' is ignored for a 2 x 2 table in: fisher.test(m, hybrid = TRUE)

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] indexing question

2006-11-13 Thread Benilton Carvalho
diff(tmp[idx])

cheers,

b

On Nov 13, 2006, at 3:06 PM, Leeds, Mark ((IED)) wrote:

 I have the following set of indices, call it idx, that correspond  
 to the
 indices of a vector say temp.

  [1]   31   36   41   61   66   71   91   96  101  121  126  131  151
 156  161  181  186  191  211  216  221  241  246  251  271  276  281
 301  306  311  331  336  341  361  366
  [36]  371  391  396  401  421  426  431  451  456  461  481  486  491
 511  516  521  541  546  551  571  576  581  601  606  611  631  636
 641  661  666  671  691  696  701  721
  [71]  726  731  751  756  761  781  786  791  811  816  821  841  846
 851  871  876  881  901  906  911  931  936  941  961  966  971  991
 996 1001 1021 1026 1031 1051 1056 1061
 [106] 1081 1086 1091  1116 1121 1141 1146 1151 1171 1176 1181 1201
 1206 1211 1231 1236 1241 1261 1266 1271 1291 1296 1301 1321 1326 1331
 1351 1356 1361 1381 1386 1391 1411 1416
 [141] 1421


 I want to calculate temp[36] - temp[31] and temp[41] - temp[36]

 Similarly, temp[66] - temp[61] and temp[71] - temp[66]
 .
 .
 .
 .
 Similarly temp[1416]-temp[1411]
   temp[1421] - temp[1416]


 I'm doing this because the above subractions represent pairs of  
 returns
 and the correlations between them wil be calculated eventually.

 In other words, eventually I will have

 X_36_31  ( i.e : temp[36] - temp[31] )
 X_66-61
 X_96-91
 .
 .
 .
 .
 .
 .
 .
 X_1411-1416

 as one vector and

 Y_41-36
 Y_71-66
 Y_101-96
 .
 .
 .
 .
 .
 Y_1416_1421

 as another vector.

 and will calculate the correlation between the two vectors in order to
 get one number.


 The point is I am  really only using the indices 31, 61, 91 etc as
 anchor's so a regular diff(temp[idx]) won't work because it will diff
 all
 the elements that are next to each other ? This is a weird problem.  
 I'm
 still thinking about it. I'm hoping to figure it out before someone
 sends me something but I won't mind so much if I get an external
 solution first. I have no pride.
 

 This is not an offer (or solicitation of an offer) to buy/se... 
 {{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
 and provide commented, minimal, self-contained, reproducible code.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] indexing question

2006-11-13 Thread Leeds, Mark \(IED\)
thanks beilton but that won't work. A diff will also include 61-41 etc
and I don't want to include those.

I'm working on using lapply or sapply with a seq along 31, 61, etc.
I'll let you know if it works.


-Original Message-
From: Benilton Carvalho [mailto:[EMAIL PROTECTED] 
Sent: Monday, November 13, 2006 3:18 PM
To: Leeds, Mark (IED)
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] indexing question

diff(tmp[idx])

cheers,

b

On Nov 13, 2006, at 3:06 PM, Leeds, Mark ((IED)) wrote:

 I have the following set of indices, call it idx, that correspond to 
 the indices of a vector say temp.

  [1]   31   36   41   61   66   71   91   96  101  121  126  131  151
 156  161  181  186  191  211  216  221  241  246  251  271  276  281
 301  306  311  331  336  341  361  366  [36]  371  391  396  401  421

 426  431  451  456  461  481  486  491
 511  516  521  541  546  551  571  576  581  601  606  611  631  636
 641  661  666  671  691  696  701  721  [71]  726  731  751  756  761

 781  786  791  811  816  821  841  846
 851  871  876  881  901  906  911  931  936  941  961  966  971  991
 996 1001 1021 1026 1031 1051 1056 1061 [106] 1081 1086 1091  1116 
 1121 1141 1146 1151 1171 1176 1181 1201
 1206 1211 1231 1236 1241 1261 1266 1271 1291 1296 1301 1321 1326 1331
 1351 1356 1361 1381 1386 1391 1411 1416 [141] 1421


 I want to calculate temp[36] - temp[31] and temp[41] - temp[36]

 Similarly, temp[66] - temp[61] and temp[71] - temp[66] .
 .
 .
 .
 Similarly temp[1416]-temp[1411]
   temp[1421] - temp[1416]


 I'm doing this because the above subractions represent pairs of  
 returns
 and the correlations between them wil be calculated eventually.

 In other words, eventually I will have

 X_36_31  ( i.e : temp[36] - temp[31] )
 X_66-61
 X_96-91
 .
 .
 .
 .
 .
 .
 .
 X_1411-1416

 as one vector and

 Y_41-36
 Y_71-66
 Y_101-96
 .
 .
 .
 .
 .
 Y_1416_1421

 as another vector.

 and will calculate the correlation between the two vectors in order to
 get one number.


 The point is I am  really only using the indices 31, 61, 91 etc as
 anchor's so a regular diff(temp[idx]) won't work because it will diff
 all
 the elements that are next to each other ? This is a weird problem.  
 I'm
 still thinking about it. I'm hoping to figure it out before someone
 sends me something but I won't mind so much if I get an external
 solution first. I have no pride.
 

 This is not an offer (or solicitation of an offer) to buy/se... 
 {{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
 and provide commented, minimal, self-contained, reproducible code.


This is not an offer (or solicitation of an offer) to buy/se...{{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
and provide commented, minimal, self-contained, reproducible code.


[R] handling time units--hours, minutes, days--from file times

2006-11-13 Thread Warren
Dear R-helpers,
I am trying to generalize my function for recording measurement times from
file times mtime--my intervals are minutes to hours over the course of
several days.  I want to use hours as my units, and I have had trouble
dealing with time units in a general way.  I have a simple solution for the
dealing with regular intervals, but I have not been able to handle irregular
intervals.  I tried using the zoo package, and had similar problems with the
conversion of time units.

Code and sample data is below--
for your info--after this, I am feeding the times and measurements to nls.

Thanks in advance,
-- Warren


#Here is the failed code:

myfilenames-row.names(myfileinfo)

fileno-length(myfilenames)

filetimes-numeric()

for (i in 2:fileno){rel.read.time-myfileinfo$mtime[i]-myfileinfo$mtime[1]

rel.read.time-(as.numeric(rel.read.time))

##this next part is problematic

##i would love to say  time.in.hours-as.hours
(myfileinfo$mtime[i]-myfileinfo$mtime[1])

ifelse(rel.read.time =5, hour.read.time-(rel.read.time /60),
hour.read.time-rel.read.time)

filetimes-c(filetimes,rel.read.time)}


#For regular intervals of readings, the code below generates measurement
times from the first interval:

myfilenames-row.names(myfileinfo)
fileno-length(myfilenames)
intervalue- myfileinfo$mtime[2] - myfileinfo$mtime[1]

unitless.interval-(as.numeric(intervalue))

ifelse(unitless.interval=5, hour.interval-(unitless.interval/60),
hour.interval-unitless.interval)

hours-0

for (i in 1:(fileno-1)){hours-c(hours,(i*hour.interval))}



myfileinfo -
size isdir mode
mtime   ctime
G2659310 2006-307-21-09-57.txt 10220 FALSE  666 2006-11-03 21:04:00
2006-11-13 10:56:39
G2659310 2006-307-21-45-55.txt 10230 FALSE  666 2006-11-03 21:40:00
2006-11-13 10:56:39
G2659310 2006-307-22-23-00.txt 10236 FALSE  666 2006-11-03 22:17:00
2006-11-13 10:56:39
G2659310 2006-307-23-00-33.txt 10236 FALSE  666 2006-11-03 22:54:00
2006-11-13 10:56:39
G2659310 2006-307-23-39-39.txt 10234 FALSE  666 2006-11-03 23:33:00
2006-11-13 10:56:39
G2659310 2006-308-00-17-25.txt 10234 FALSE  666 2006-11-04 00:11:00
2006-11-13 10:56:39
G2659310 2006-308-00-53-45.txt 10228 FALSE  666 2006-11-04 00:48:00
2006-11-13 10:56:39
G2659310 2006-308-01-30-03.txt 10208 FALSE  666 2006-11-04 01:25:00
2006-11-13 10:56:39
G2659310 2006-308-02-06-00.txt 10188 FALSE  666 2006-11-04 02:00:00
2006-11-13 10:56:39
G2659310 2006-308-02-41-57.txt 10168 FALSE  666 2006-11-04 02:36:00
2006-11-13 10:56:39
G2659310 2006-308-03-17-54.txt 10158 FALSE  666 2006-11-04 03:12:00
2006-11-13 10:56:39
G2659310 2006-308-03-53-51.txt 10148 FALSE  666 2006-11-04 03:48:00
2006-11-13 10:56:39
G2659310 2006-308-04-29-48.txt 10144 FALSE  666 2006-11-04 04:24:00
2006-11-13 10:56:39
G2659310 2006-308-05-05-44.txt 10132 FALSE  666 2006-11-04 05:00:00
2006-11-13 10:56:39
G2659310 2006-308-05-41-41.txt 10136 FALSE  666 2006-11-04 05:35:00
2006-11-13 10:56:39
G2659310 2006-308-06-17-39.txt 10128 FALSE  666 2006-11-04 06:11:00
2006-11-13 10:56:39
G2659310 2006-308-06-52-27.txt 10128 FALSE  666 2006-11-04 06:47:00
2006-11-13 10:56:39
G2659310 2006-308-07-28-25.txt 10130 FALSE  666 2006-11-04 07:23:00
2006-11-13 10:56:39
G2659310 2006-308-08-04-22.txt 10130 FALSE  666 2006-11-04 07:59:00
2006-11-13 10:56:39
G2659310 2006-308-08-40-20.txt 10134 FALSE  666 2006-11-04 08:35:00
2006-11-13 10:56:39
G2659310 2006-308-09-16-18.txt 10166 FALSE  666 2006-11-04 09:11:00
2006-11-13 10:56:39
G2659310 2006-308-09-52-15.txt 10180 FALSE  666 2006-11-04 09:46:00
2006-11-13 10:56:39
G2659310 2006-308-10-28-13.txt 10192 FALSE  666 2006-11-04 10:22:00
2006-11-13 10:56:39
G2659310 2006-308-11-04-10.txt 10210 FALSE  666 2006-11-04 10:58:00
2006-11-13 10:56:39
G2659310 2006-308-11-40-08.txt 10210 FALSE  666 2006-11-04 11:34:00
2006-11-13 10:56:39
G2659310 2006-308-12-16-06.txt 10222 FALSE  666 2006-11-04 12:10:00
2006-11-13 10:56:39
G2659310 2006-308-12-50-55.txt 10222 FALSE  666 2006-11-04 12:46:00
2006-11-13 10:56:39
G2659310 2006-308-13-26-54.txt 10224 FALSE  666 2006-11-04 13:21:00
2006-11-13 10:56:39
G2659310 2006-308-14-02-52.txt 10220 FALSE  666 2006-11-04 13:57:00
2006-11-13 10:56:39
G2659310 2006-308-14-38-51.txt 10224 FALSE  666 2006-11-04 14:33:00
2006-11-13 10:56:39
G2659310 2006-308-15-14-50.txt 10216 FALSE  666 2006-11-04 15:09:00
2006-11-13 10:56:39
G2659310 2006-308-15-50-47.txt 10202 FALSE  666 2006-11-04 15:45:00
2006-11-13 10:56:39
G2659310 2006-308-16-26-45.txt 10202 FALSE  666 2006-11-04 16:21:00
2006-11-13 10:56:39
G2659310 2006-308-17-02-43.txt 10190 FALSE  666 2006-11-04 16:57:00
2006-11-13 10:56:39
G2659310 2006-308-17-38-41.txt 10188 FALSE  666 2006-11-04 17:32:00
2006-11-13 10:56:39
G2659310 2006-308-18-14-40.txt 10174 FALSE  666 2006-11-04 18:08:00
2006-11-13 10:56:39
G2659310 2006-308-18-49-29.txt 10170 FALSE  666 2006-11-04 18:44:00
2006-11-13 10:56:39
G2659310 2006-308-19-25-28.txt 10158 FALSE  666 2006-11-04 19:20:00
2006-11-13 10:56:39
G2659310 

Re: [R] indexing question

2006-11-13 Thread Ray Brownrigg
On Tuesday 14 November 2006 09:28, Leeds, Mark (IED) wrote:
 thanks beilton but that won't work. A diff will also include 61-41 etc
 and I don't want to include those.

 I'm working on using lapply or sapply with a seq along 31, 61, etc.
 I'll let you know if it works.


Try looking at:
dim(idx) - c(3, length(idx)/3)
then:
tmp[idx[2, ]] - temp[idx[1, ]]
temp[idx[3, ]] - temp[idx[2, ]]

HTH
Ray Brownrigg

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] R and Fortran 9x -- advice

2006-11-13 Thread Mike Prager
Christos Hatzis [EMAIL PROTECTED] wrote:

 Can you recommend any good books on Fortran 90/95?  I had been an old user
 of Fortran 77 but haven't followed the developments in the last 15 years or
 so...
 
Ravi Varadhan has already recommended Metcalf and Reid -
Fortran 90/95 Explained.  I own and like that book, so I'll
second that, with a few caveats--

1. The same authors + M. Cohen have released Fortran 95/2003
Explained.  If you buy, why not buy the latest?  There are as
yet no Fortran 2003 compilers, but F95 compilers are starting to
introduce F2003 features.  

2. The style of MR(C) is quite condensed.  Think R help pages.
If you respond to terse explanations, this is your book. If not,
check Amazon.com.

3. The book is NOT specifically written to help those moving
from F77 to F95.  For that, you might look at a Web tutorial,
such as

http://www.nsc.liu.se/~boein/f77to90/f77to90.html


--and a couple of comments:

Because standard F95, unlike F77, meets most programming needs,
many F95 programmers try to write strictly standard-conforming
code, rather than use vendor extensions.  (I have ported a
10,000+ line program from Windows to Linux with NO changes.) A
book with a detailed explanation of the Fortran 95 standard is
Fortran 95 Handbook by Adams et al.  I have a copy and use it
rarely, but when I just can't seem to find something in MRC,
I'm glad to have it.

The Usenet group comp.lang.fortran has a high S/N ratio and
courteous tone, with several members of the Fortran standards
committee being frequent contributors.  It's a great resource.

HTH,
...Mike

-- 
Mike Prager, NOAA, Beaufort, NC
* Opinions expressed are personal and not represented otherwise.
* Any use of tradenames does not constitute a NOAA endorsement.

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Nominal Respose Model in R

2006-11-13 Thread Xavier Giovanni Ordóñez Camacho
Hi:

I have been working in Item Response Theory, exactly, with Nominal Response 
Model (NRM). Exist in R 
a function for estimate parameter and ability from database for this Model?.

Thank you,


Xavier G. Ordóñez



[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] question on MSM warning message

2006-11-13 Thread Im, Kelly
Hello

After (simple random cluster) resampling with replacement I ran MSM
function and I'm getting the following warning message ,which I'm not
sure why. I don't have any absorbing stage set up in my MSM model. I
have a 4 stage uni-directional MSM model. The only thing I see might be
a problem is that when the last stage (stage 4 in my data) gets repeated
it seems to give me a warning message. although, it picks only one
observation whereas there might be other such cases in the data
(printing only the first case?) 

 out1.msm2
-msm(stage~years,subject=patient,data=mydata2,qmatrix=new.twoway9.q,con
trol=list(trace=1,report=1),covariates=~black+tport11,constraint=list(bl
ack=c(1,1,1),tport11=c(1,1,1)))
...
Warning message:
Absorbing - absorbing transition at observation 30 in:
msm.check.model(msmdata$fromstate, msmdata$tostate, msmdata$obs,  

My data at observation around 30 look like ('years' is the time
variable)

  mydata2[27:35,]
count patient black years stage tport11 pcount
 11029 1 0 1   1 17
 1   1029 1 0 1   1 17
 2   1029 1194   1 17
 2   1029 1194   1 17
 1   1033 0 0 1   1 20
 1   1033 0 0 1   1 20
 2   1033 0184   1 20
 2   1033 0184   1 20

I'm wondering if I can ignore this warning message or should I have not
done resampling with replacement. Any input would be greatly
appreciated.

Thank you so much for your time,

kelly
-
KyungAh (Kelly)  Im   
Epidemiology Data Center 
Graduate School of Public Health
Department of Epidemiology
University of Pittsburgh
PA 15261
PHONE:(412) 624- 4612 FAX: (412) 624- 3775
Email: [EMAIL PROTECTED]
http://www.edc.gsph.pitt.edu

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] For MacBook, best way to do R?

2006-11-13 Thread Mitchell Maltenfort
I'm getting a MacBook and I'd like to stick with OS X rather than
convert it to Linux just yet.

However, my main concern is having decent performance.

What's my best option:

*use the existing binary for R?

*compile R fresh under OS X?

* install Linux and run R under that?

Does anyone have any recent experience they can share?  Thanks!

-- 
I can answer any question.
I don't know is an answer.
I don't know yet is a better answer.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] A printing macro

2006-11-13 Thread Murray Jorgensen

Thanks for these suggestions, Professor Ripley. It's interesting that
the function parameters in R are not truly dummy as they can effect
the result of a function.

Murray

Prof Brian Ripley wrote:
 ddisp - function(dvar) {
 yn - substitute(dvar)
 csqt - eval.parent(substitute(chisq.test(G,dvar), list(dvar=yn)))
 
 }
 
 There are other ways, such as forming the cross-classification table, 
 setting its dimnames and passing that to chisq.test.
 
 On Mon, 13 Nov 2006, Murray Jorgensen wrote:
 
 I am exploring the result of clustering a large multivariate data set
 into a number of groups, represented, say, by a factor G.

 I wrote a function to see how categorical variables vary between groups:

   ddisp - function(dvar) {
 +  csqt - chisq.test(G,dvar)
 +  print(csqt$statistic)
 +  print(csqt$observed)
 +  print(round(csqt$expected))
 +  round(csqt$residuals)
 +  }
 
   x - ceiling(4*runif(100))
   G - gl(4,1,100)
   ddisp(x)
 X-squared
  6.235645
dvar
 G1  2  3  4
   1 10  5  5  5
   2  6  9  5  5
   3  8  6  5  6
   4  7  4  4 10
dvar
 G   1 2 3 4
   1 8 6 5 6
   2 8 6 5 6
   3 8 6 5 6
   4 8 6 5 6
dvar
 G1  2  3  4
   1  1  0  0 -1
   2 -1  1  0 -1
   3  0  0  0  0
   4  0 -1  0  1
 Warning message:
 Chi-squared approximation may be incorrect in: chisq.test(G, dvar)

 As I need to apply this function to a large number of variables x it
 would be helpful if the function printed x rather than the formal
 argument dvar. I have a vague idea that things like deparse() and
 substitute() will come into the solution but I have not yet come up with
 the right incantation. Any help appreciated!

 Murray Jorgensen


 

-- 
Dr Murray Jorgensen  http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: [EMAIL PROTECTED]Fax 7 838 4155
Phone  +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862


-- 
Dr Murray Jorgensen  http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: [EMAIL PROTECTED]Fax 7 838 4155
Phone  +64 7 838 4773 wkHome +64 7 825 0441Mobile 021 1395 862

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Converting Character Strings to Numerical Values

2006-11-13 Thread Peter Lauren
Does R have anything like strtod, atoi, atod, etc. in
C?  I would like to covert 3 to 3 and 3.5 to 3.5 .

Thanks,
Peter Lauren.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Nominal Respose Model in R

2006-11-13 Thread Dimitrios Rizopoulos
AFAIK there exist functions to fit IRT models only for ordinal data.

In particular, look at function grm() that fits the Graded Response  
Model in the 'ltm' package (using Marginal Maximum Likelihood),  
function PCM() in the 'eRm' package (using Condtional Maximum  
Likelihood), and function MCMCordfactanal() in the 'MCMCpack' package  
(using a Bayesian approach).

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/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
  http://www.student.kuleuven.be/~m0390867/dimitris.htm


Quoting Xavier Giovanni Ordóñez Camacho [EMAIL PROTECTED]:

 Hi:

 I have been working in Item Response Theory, exactly, with Nominal   
 Response Model (NRM). Exist in R
 a function for estimate parameter and ability from database for this Model?.

 Thank you,


 Xavier G. Ordóñez



   [[alternative HTML version deleted]]





Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm

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


Re: [R] For MacBook, best way to do R?

2006-11-13 Thread Lanre Okusanya
Running R using the precomplied binary works really well and you
should have no issues. Unless you really know what you are doing
compiling R to work on a mac from sctrach can be at the very least
tedious, at most, daunting, and there is really no point in doing so,
since the binary already takes care of all the issues you may run
into. In addition, you will need in install additional software (free
though, from the OS X dvd) before you can compile.


Lanre

On 11/13/06, Mitchell Maltenfort [EMAIL PROTECTED] wrote:
 I'm getting a MacBook and I'd like to stick with OS X rather than
 convert it to Linux just yet.

 However, my main concern is having decent performance.

 What's my best option:

 *use the existing binary for R?

 *compile R fresh under OS X?

 * install Linux and run R under that?

 Does anyone have any recent experience they can share?  Thanks!

 --
 I can answer any question.
 I don't know is an answer.
 I don't know yet is a better answer.

 __
 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
 and provide commented, minimal, self-contained, reproducible code.


__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Converting Character Strings to Numerical Values

2006-11-13 Thread Sarah Goslee
?as.numeric

 x - 3
 y - 3.5

 x
[1] 3
 as.numeric(x)
[1] 3

 y
[1] 3.5
 as.numeric(y)
[1] 3.5

Sarah

On 11/13/06, Peter Lauren [EMAIL PROTECTED] wrote:
 Does R have anything like strtod, atoi, atod, etc. in
 C?  I would like to covert 3 to 3 and 3.5 to 3.5 .

 Thanks,
 Peter Lauren.

 __
 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
 and provide commented, minimal, self-contained, reproducible code.



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

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Multivariate time-series

2006-11-13 Thread David Kaplan
Hi all,

I'm looking for R packages that estimate multivariate time-series models 
or vector-autoregression (VAR) time-series models.

Thanks

David


-- 
===
David Kaplan, Ph.D.
Professor
Department of Educational Psychology
University of Wisconsin - Madison
Educational Sciences, Room, 1061
1025 W. Johnson Street
Madison, WI 53706

email: [EMAIL PROTECTED]
homepage: http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm
Phone: 608-262-0836

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Multivariate time-series

2006-11-13 Thread Leeds, Mark \(IED\)
I think I remember seeing something called sem ? It's listed as on the
packages on
www.r-project.org and the explanation is what you want. I just can't be
sure of the name.
There are probably more than just sem.

  

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of David Kaplan
Sent: Monday, November 13, 2006 4:42 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Multivariate time-series

Hi all,

I'm looking for R packages that estimate multivariate time-series models
or vector-autoregression (VAR) time-series models.

Thanks

David


--

===
David Kaplan, Ph.D.
Professor
Department of Educational Psychology
University of Wisconsin - Madison
Educational Sciences, Room, 1061
1025 W. Johnson Street
Madison, WI 53706

email: [EMAIL PROTECTED]
homepage:
http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm
Phone: 608-262-0836

__
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
and provide commented, minimal, self-contained, reproducible code.


This is not an offer (or solicitation of an offer) to buy/se...{{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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] hybrid in fisher.test broken?

2006-11-13 Thread Peter Dalgaard
Jonathan Dushoff [EMAIL PROTECTED] writes:

 The hybrid feature in fisher.test looks to me like an excellent way to
 analyze my two-way tables.   The only problem is that it does not seem
 to be implemented.  Am I right about this?
 
 An example is pasted below.  I note that I get the warning message only
 when I shouldn't: for a 2x2 table hybrid seems to be ignored without
 warning.  In no case does fisher.test seem to be checking Cochran
 criteria as promised.
 
 Any help will be appreciated.


Well, something seems to have been mangled in connection with these
changes: 


r36358 | ripley | 2005-11-15 17:50:12 +0100 (Tue, 15 Nov 2005) | 2
lines

add simulation option to fisher.test


r36266 | ripley | 2005-11-10 22:19:34 +0100 (Thu, 10 Nov 2005) | 2
lines

move 2x2 case in fisher.test to similar code as that used for or != 1



so that the warning 

else if (hybrid) {
warning('hybrid' is ignored for a 2 x 2 table)

now appears inside

if (nr != 2 || nc != 2) {

which is clearly wrong.

However, the call to fexact that follows does have parameters set
for the hybrid algorithm, so it's just the warning that is spurious.
Notice that the hybrid algorithm may save a bit of time on probability
calculations, but still needs to generate a large number of tables, so
you really do need a larger workspace:

 fisher.test(matrix(10*1:9, nc=3),hybrid=TRUE,workspace=1e6)

Fisher's Exact Test for Count Data

data:  matrix(10 * 1:9, nc = 3) 
p-value = 0.3149
alternative hypothesis: two.sided 

Warning message:
'hybrid' is ignored for a 2 x 2 table in: fisher.test(matrix(10 * 1:9,
nc = 3), hybrid = TRUE, workspace = 1e+06) 

whereas (notice the slightly different p value)

 fisher.test(matrix(10*1:9, nc=3),hybrid=FALSE,workspace=1e6)

Fisher's Exact Test for Count Data

data:  matrix(10 * 1:9, nc = 3) 
p-value = 0.3145
alternative hypothesis: two.sided 






 JD
 
 --
 
 
 R : Copyright 2006, The R Foundation for Statistical Computing
 Version 2.3.1 (2006-06-01)
 
 .
 .
 .
 
  m = matrix(10*1:9, nc=3)
  fisher.test(m, hybrid=TRUE)
 Error in fisher.test(m, hybrid = TRUE) : FEXACT error 6.
 LDKEY is too small for this problem.
 Try increasing the size of the workspace.
 In addition: Warning message:
 'hybrid' is ignored for a 2 x 2 table in: fisher.test(m, hybrid = TRUE)
 
 __
 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
 and provide commented, minimal, self-contained, reproducible code.
 

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

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


[R] Creating data for logistic regression and Cox proportional hazards regression

2006-11-13 Thread John Sorkin
I know that mvrnorm from MASS (generously provided by Profs. Venables
and Ripley) can be used to generate multivariable normal data that can
be used in a linear regression with certain desired characteristics
(e.g. a given mean for each variable as well as a given
variance-covariance pattern). Is there any similar facility that can be
used to generate data for (1) a logistic regression and (2) a Cox
proportional hazards regression?  
Thanks,
John
 
 
John Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
Baltimore VA Medical Center GRECC,
University of Maryland School of Medicine Claude D. Pepper OAIC,
University of Maryland Clinical Nutrition Research Unit, and
Baltimore VA Center Stroke of Excellence

University of Maryland School of Medicine
Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524

(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
[EMAIL PROTECTED]
Confidentiality Statement:
This email message, including any attachments, is for the so...{{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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Multivariate time-series

2006-11-13 Thread Achim Zeileis
On Mon, 13 Nov 2006 15:42:06 -0600 David Kaplan wrote:

 Hi all,
 
 I'm looking for R packages that estimate multivariate time-series
 models or vector-autoregression (VAR) time-series models.

The Econometrics task view at
  http://CRAN.R-project.org/src/contrib/Views/Econometrics.html
has in the Time-series modelling section:

For estimating VAR models, several methods are available: simple models
can be fitted by ar() in stats, more elaborate models are provided
package vars, estVARXls() in dse and a Bayesian approach is
available in MSBVAR.

hth,
Z

 Thanks
 
 David
 
 
 -- 
 ===
 David Kaplan, Ph.D.
 Professor
 Department of Educational Psychology
 University of Wisconsin - Madison
 Educational Sciences, Room, 1061
 1025 W. Johnson Street
 Madison, WI 53706
 
 email: [EMAIL PROTECTED]
 homepage:
 http://www.education.wisc.edu/edpsych/facstaff/kaplan/kaplan.htm
 Phone: 608-262-0836
 
 __
 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 and provide commented,
 minimal, self-contained, reproducible code.


__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Embedded carriage returns in text document

2006-11-13 Thread Dennis Fisher
Colleagues,

I am using R 2.4.0 on both a Mac (10.4.8) and Linux (RedHat 9).  To  
read data from an Excel spreadsheet, I do save as  in Excel, then  
select the Text (tab-delimited) format.  The resulting file uses a  
tab separator and I can usually read the file using read.delim.

Sometimes, the header row contains embedded carriage returns.  When I  
view the file, these carriage returns appear as  ^M.

Now the problem:
When I read.delim these files, they do not read correctly.  Sometimes  
I get error messages; sometimes only the first line is read.   
Interestingly, invoking the option skip=1 (or a larger N) does not  
appear to bypass the problem.

I can solve the problem by manually deleting these carriage returns  
either in the original Excel file or the .txt version.   However,  
this is not an ideal solution.

Does anyone have a work-around within R?

Dennis


Dennis Fisher MD
P  (The P Less Than Company)
Phone: 1-866-PLessThan (1-866-753-7784)
Fax: 1-415-564-2220
www.PLessThan.com   


[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Creating data for logistic regression and Cox proportional hazards regression

2006-11-13 Thread Frank E Harrell Jr
John Sorkin wrote:
 I know that mvrnorm from MASS (generously provided by Profs. Venables
 and Ripley) can be used to generate multivariable normal data that can
 be used in a linear regression with certain desired characteristics
 (e.g. a given mean for each variable as well as a given
 variance-covariance pattern). Is there any similar facility that can be
 used to generate data for (1) a logistic regression and (2) a Cox
 proportional hazards regression?  
 Thanks,
 John

See the examples in the help file for the lrm and cph functions in the 
Design package as well as in the spower function in the Hmisc package.

Frank

  
  
 John Sorkin M.D., Ph.D.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] hybrid in fisher.test broken?

2006-11-13 Thread Jonathan Dushoff
Thank you very much for this useful response.  

I did a search for fexact both inside R and on R-project.org before
posting: can somebody tell me where I could have found this information?

In particular, I am still not clear on whether the Cochran criteria are
being tested.

Thanks,

JD

--

Jonathan Dushoff [EMAIL PROTECTED] writes:

 The hybrid feature in fisher.test looks to me like an excellent way to
 analyze my two-way tables.   The only problem is that it does not seem
 to be implemented.  Am I right about this?
 
 An example is pasted below.  I note that I get the warning message only
 when I shouldn't: for a 2x2 table hybrid seems to be ignored without
 warning.  In no case does fisher.test seem to be checking Cochran
 criteria as promised.
 
 Any help will be appreciated.


Well, something seems to have been mangled in connection with these
changes: 


r36358 | ripley | 2005-11-15 17:50:12 +0100 (Tue, 15 Nov 2005) | 2
lines

add simulation option to fisher.test


r36266 | ripley | 2005-11-10 22:19:34 +0100 (Thu, 10 Nov 2005) | 2
lines

move 2x2 case in fisher.test to similar code as that used for or != 1



so that the warning 

else if (hybrid) {
warning('hybrid' is ignored for a 2 x 2 table)

now appears inside

if (nr != 2 || nc != 2) {

which is clearly wrong.

However, the call to fexact that follows does have parameters set
for the hybrid algorithm, so it's just the warning that is spurious.
Notice that the hybrid algorithm may save a bit of time on probability
calculations, but still needs to generate a large number of tables, so
you really do need a larger workspace:

 fisher.test(matrix(10*1:9, nc=3),hybrid=TRUE,workspace=1e6)

Fisher's Exact Test for Count Data

data:  matrix(10 * 1:9, nc = 3) 
p-value = 0.3149
alternative hypothesis: two.sided 

Warning message:
'hybrid' is ignored for a 2 x 2 table in: fisher.test(matrix(10 * 1:9,
nc = 3), hybrid = TRUE, workspace = 1e+06) 

whereas (notice the slightly different p value)

.
.
.

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Fitting mean and covariance of Multivariate normal with censored data

2006-11-13 Thread Sicotte, Hugues Ph.D.
Nobody answered, but this is what I did.

I used an iterative poor man's imputation and filled in the data with a
single value at each iteration. My covariance matrices and variance will
be a little underestimated, but it'll do for me at this stage of this
project.

I iteratively filled in the data using the current data to estimate the
parameters of the normal distribution.
Using the current estimates, I figured out how much of the current
distribution is outside each cut-off. Then fill-in at the mid-point of
the probability density that is outside the area. I limited my
corrections to 3*sigma.

Here is the code section inside the loop (I am not showing another part
of the code that does outlier elimination, nor the rest of the loop that
recomputes the covariance and the parameters of the normal distribution)

Trim the data to only use lines with at least 60% real data.
Trim the columns to those that have enough data to estimate the
parameters of the normal distribution (e.g. ~ 20)

Loop over iter iterations ( from 1 to iterlast)

  Compute mean, variance, and co-variance matrix.
  Code to detect outliers.. Only activated after
iterlast iterations
   Sort and trim worst outliers (left and
right, up to 10% of data)

  For Each line of data, find out which elements were
missing in original data set.

for(im in missings) {
   # Compute the quantile on the cutoff away
from N(mode,variance) for both cutoff
 
quantOver-pnorm(c(highR),mean=mu[im],sd=sigmaMu[im])
 
quantUnder-pnorm(c(lowR),mean=mu[im],sd=sigmaMu[im])
   # Find out which end of the distribution is
the most outside of the limit
   quantOverMid- (1.0-quantOver)/2
   quantUnderMid- quantUnder/2
#  correct the mean, but underestimate the
variance a little bit in cases
#  where width of pattern width of
cutoff-region
   fixUp-mu[im]
   fixLow-mu[im]
#  Don't fix if distribution is cutoff at edges.
   if(quantOverMid0.05) {
quantOverMid-0.0
   } else {
 
fixUp-qnorm(c(1.0-quantOverMid),mean=mu[im],sd=sigmaMu[im])
   }
   if(quantUnderMid0.05) {
quantUnderMid-0.0
   } else {
 
fixLow-qnorm(c(quantUnderMid),mean=mu[im],sd=sigmaMu[im])
   }
#  Censored variables pull from both sides.
   quantNormSum - (quantOverMid+quantUnderMid)
   if(quantNormSum=0.05) {
   quantNormSum- 1.0/quantNormSum
   fix -
(fixUp*quantOverMid+fixLow*quantUnderMid)*quantNormSum
   } else {
#  If missing data is in the tails, correct
by cutoff.
   if(quantNormSum1e-10) {
   fix-
(highR*(1.0-quantOver)+lowR*quantUnder)/((1.0-quantOver)+quantUnder)
   } else {
   fix-mu[im]
   }
   }
   if(!is.finite(fix)) {
   if(1.0-quantOver  quantUnder) {
fix-mu[im]+3*sigmaMu[im]
   } else {
 fix-  mu[im]-3*sigmaMu[im]
   }
   } else {
   if(fixmu[im]+3*sigmaMu[im]) {
 fix-mu[im]+3*sigmaMu[im]
   } else if (fix mu[im]-3*sigmaMu[im]) {
 fix-  mu[im]-3*sigmaMu[im]
   }
   }




 _ 
 From: Sicotte, Hugues   Ph.D.  
 Sent: Wednesday, November 01, 2006 11:38 AM
 To:   'r-help@stat.math.ethz.ch'
 Subject:  Fitting mean and covariance of Multivariate normal with
 censored data
 
 Hello,
 I have been googling for 2 days and I cannot find the answer
 in previous posts.
 
 I have a set of d-dimensional data elements (d=11 .. 14), each data
 point can be censored at different values both
 Lower-limit and upper limit.
 
 N = 2000 sets of vectors of
 D=11 data points per vector.
 Each of the N*D points can have different upper and lower limits.
 
 I simply want to fit a Multivariate distribution to the data and get
 the mean and covariance matrix  of the sample.
 
 The closest post I saw mentionned the function survreg in the
 survival5 package, but I can't figure out how
 To use it for a 

Re: [R] Embedded carriage returns in text document

2006-11-13 Thread Peter Dalgaard
Dennis Fisher [EMAIL PROTECTED] writes:

 Colleagues,
 
 I am using R 2.4.0 on both a Mac (10.4.8) and Linux (RedHat 9).  To  
 read data from an Excel spreadsheet, I do save as  in Excel, then  
 select the Text (tab-delimited) format.  The resulting file uses a  
 tab separator and I can usually read the file using read.delim.
 
 Sometimes, the header row contains embedded carriage returns.  When I  
 view the file, these carriage returns appear as  ^M.
 
 Now the problem:
 When I read.delim these files, they do not read correctly.  Sometimes  
 I get error messages; sometimes only the first line is read.   
 Interestingly, invoking the option skip=1 (or a larger N) does not  
 appear to bypass the problem.
 
 I can solve the problem by manually deleting these carriage returns  
 either in the original Excel file or the .txt version.   However,  
 this is not an ideal solution.
 
 Does anyone have a work-around within R?

Hmm,... I suppose that CR messes with what R or the system thinks is
the line-end character in this particular file.

My first idea would be to read from a pipe() which executed 

sed 's/\r//' myfile.dat

or something in that vein. Beware of quoting and differences in sed
versions  



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

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


Re: [R] handling time units--hours, minutes, days--from file times

2006-11-13 Thread Tim Calkins
hardly the most efficient way to go, but consider using a substring function
to extract the time bits from your data, then reading them as POSIX dates
and using difftime.

for example,

 mytime - c(G2659310 2006-310-10-55-32.txt 10134 FALSE  666 2006-11-06
10:49:00 2006-11-13 10:56:41,
G2659310 2006-310-11-31-43.txt 10136 FALSE  666 2006-11-06 11:25:00
2006-11-13 10:56:41,
G2659310 2006-310-12-08-04.txt 10140 FALSE  666 2006-11-06 12:02:00
2006-11-13 10:56:41,
G2659310 2006-310-12-43-17.txt 10148 FALSE  666 2006-11-06 12:38:00
2006-11-13 10:56:41,
G2659310 2006-310-13-20-21.txt 10138 FALSE  666 2006-11-06 13:15:00
2006-11-13 10:56:41,
G2659310 2006-310-13-56-28.txt 10080 FALSE  666 2006-11-06 13:51:00
2006-11-13 10:56:38)
 time1 - substring(mytime,49,67)
 time2 - substring(mytime,69,87)

interval - difftime(as.POSIXct(time1, tz = AET10ADT), as.POSIXct(time2,
tz = AET10ADT), units = hours)
interval
   Time differences of -168.1281, -167.5281, -166.9114, -166.3114, -165.6947,
-165.0939 hours

i'm not exactly sure what you're after, but here's a vector of time
differences.  in general the POSIX class holds the most detailed information
about time -- note the specification of tz (time zone) -- and should be used
when a great deal of precision is required.

tc



On 11/14/06, Warren [EMAIL PROTECTED] wrote:

 Dear R-helpers,
 I am trying to generalize my function for recording measurement times from
 file times mtime--my intervals are minutes to hours over the course of
 several days.  I want to use hours as my units, and I have had trouble
 dealing with time units in a general way.  I have a simple solution for
 the
 dealing with regular intervals, but I have not been able to handle
 irregular
 intervals.  I tried using the zoo package, and had similar problems with
 the
 conversion of time units.

 Code and sample data is below--
 for your info--after this, I am feeding the times and measurements to nls.

 Thanks in advance,
 -- Warren


 #Here is the failed code:

 myfilenames-row.names(myfileinfo)

 fileno-length(myfilenames)

 filetimes-numeric()

 for (i in 2:fileno){rel.read.time-myfileinfo$mtime[i]-myfileinfo$mtime[1]

 rel.read.time-(as.numeric(rel.read.time))

 ##this next part is problematic

 ##i would love to say  time.in.hours-as.hours
 (myfileinfo$mtime[i]-myfileinfo$mtime[1])

 ifelse(rel.read.time =5, hour.read.time-(rel.read.time /60),
 hour.read.time-rel.read.time)

 filetimes-c(filetimes,rel.read.time)}


 #For regular intervals of readings, the code below generates measurement
 times from the first interval:

 myfilenames-row.names(myfileinfo)
 fileno-length(myfilenames)
 intervalue- myfileinfo$mtime[2] - myfileinfo$mtime[1]

 unitless.interval-(as.numeric(intervalue))

 ifelse(unitless.interval=5, hour.interval-(unitless.interval/60),
 hour.interval-unitless.interval)

 hours-0

 for (i in 1:(fileno-1)){hours-c(hours,(i*hour.interval))}



 myfileinfo -
 size isdir mode
 mtime   ctime
 G2659310 2006-307-21-09-57.txt 10220 FALSE  666 2006-11-03 21:04:00
 2006-11-13 10:56:39
 G2659310 2006-307-21-45-55.txt 10230 FALSE  666 2006-11-03 21:40:00
 2006-11-13 10:56:39
 G2659310 2006-307-22-23-00.txt 10236 FALSE  666 2006-11-03 22:17:00
 2006-11-13 10:56:39
 G2659310 2006-307-23-00-33.txt 10236 FALSE  666 2006-11-03 22:54:00
 2006-11-13 10:56:39
 G2659310 2006-307-23-39-39.txt 10234 FALSE  666 2006-11-03 23:33:00
 2006-11-13 10:56:39
 G2659310 2006-308-00-17-25.txt 10234 FALSE  666 2006-11-04 00:11:00
 2006-11-13 10:56:39
 G2659310 2006-308-00-53-45.txt 10228 FALSE  666 2006-11-04 00:48:00
 2006-11-13 10:56:39
 G2659310 2006-308-01-30-03.txt 10208 FALSE  666 2006-11-04 01:25:00
 2006-11-13 10:56:39
 G2659310 2006-308-02-06-00.txt 10188 FALSE  666 2006-11-04 02:00:00
 2006-11-13 10:56:39
 G2659310 2006-308-02-41-57.txt 10168 FALSE  666 2006-11-04 02:36:00
 2006-11-13 10:56:39
 G2659310 2006-308-03-17-54.txt 10158 FALSE  666 2006-11-04 03:12:00
 2006-11-13 10:56:39
 G2659310 2006-308-03-53-51.txt 10148 FALSE  666 2006-11-04 03:48:00
 2006-11-13 10:56:39
 G2659310 2006-308-04-29-48.txt 10144 FALSE  666 2006-11-04 04:24:00
 2006-11-13 10:56:39
 G2659310 2006-308-05-05-44.txt 10132 FALSE  666 2006-11-04 05:00:00
 2006-11-13 10:56:39
 G2659310 2006-308-05-41-41.txt 10136 FALSE  666 2006-11-04 05:35:00
 2006-11-13 10:56:39
 G2659310 2006-308-06-17-39.txt 10128 FALSE  666 2006-11-04 06:11:00
 2006-11-13 10:56:39
 G2659310 2006-308-06-52-27.txt 10128 FALSE  666 2006-11-04 06:47:00
 2006-11-13 10:56:39
 G2659310 2006-308-07-28-25.txt 10130 FALSE  666 2006-11-04 07:23:00
 2006-11-13 10:56:39
 G2659310 2006-308-08-04-22.txt 10130 FALSE  666 2006-11-04 07:59:00
 2006-11-13 10:56:39
 G2659310 2006-308-08-40-20.txt 10134 FALSE  666 2006-11-04 08:35:00
 2006-11-13 10:56:39
 G2659310 2006-308-09-16-18.txt 10166 FALSE  666 2006-11-04 09:11:00
 2006-11-13 10:56:39
 G2659310 2006-308-09-52-15.txt 10180 FALSE  666 2006-11-04 09:46:00
 2006-11-13 10:56:39
 G2659310 

Re: [R] hybrid in fisher.test broken?

2006-11-13 Thread Peter Dalgaard
Jonathan Dushoff [EMAIL PROTECTED] writes:

 Thank you very much for this useful response.  
 
 I did a search for fexact both inside R and on R-project.org before
 posting: can somebody tell me where I could have found this information?
 
 In particular, I am still not clear on whether the Cochran criteria are
 being tested.

You need the actual source for FEXACT for this, then look at the 5th
argument:

EXPECT  - Expected value used in the hybrid algorithm for
  deciding when to use asymptotic theory
  probabilities.
  (Input)
  If EXPECT = 0.0 then asymptotic theory probabilities
  are not used and Fisher exact test probabilities are
  computed.  Otherwise, if PERCNT or more of the cells in
  the remaining table have estimated expected values of
  EXPECT or more, with no remaining cell having expected
  value less than EMIN, then asymptotic chi-squared
  probabilities are used.  See the algorithm section of
  the
  manual document for details.
  Use EXPECT = 5.0 to obtain the 'Cochran' condition.

and compare with

PVAL - .C(fexact, nr, nc, x, nr, as.double(5), 
as.double(80), as.double(1), double(1), p = double(1), 
as.integer(workspace), mult = as.integer(mult), 
PACKAGE = stats)$p

The source file in question is

 src/library/stats/src/fexact.c 

available in the R source tarballs or via https://svn.r-project.org/R

 
 Thanks,
 
 JD
 
 --
 
 Jonathan Dushoff [EMAIL PROTECTED] writes:
 
  The hybrid feature in fisher.test looks to me like an excellent way to
  analyze my two-way tables.   The only problem is that it does not seem
  to be implemented.  Am I right about this?
  
  An example is pasted below.  I note that I get the warning message only
  when I shouldn't: for a 2x2 table hybrid seems to be ignored without
  warning.  In no case does fisher.test seem to be checking Cochran
  criteria as promised.
  
  Any help will be appreciated.
 
 
 Well, something seems to have been mangled in connection with these
 changes: 
 
 
 r36358 | ripley | 2005-11-15 17:50:12 +0100 (Tue, 15 Nov 2005) | 2
 lines
 
 add simulation option to fisher.test
 
 
 r36266 | ripley | 2005-11-10 22:19:34 +0100 (Thu, 10 Nov 2005) | 2
 lines
 
 move 2x2 case in fisher.test to similar code as that used for or != 1
 
 
 
 so that the warning 
 
 else if (hybrid) {
 warning('hybrid' is ignored for a 2 x 2 table)
 
 now appears inside
 
 if (nr != 2 || nc != 2) {
 
 which is clearly wrong.
 
 However, the call to fexact that follows does have parameters set
 for the hybrid algorithm, so it's just the warning that is spurious.
 Notice that the hybrid algorithm may save a bit of time on probability
 calculations, but still needs to generate a large number of tables, so
 you really do need a larger workspace:
 
  fisher.test(matrix(10*1:9, nc=3),hybrid=TRUE,workspace=1e6)
 
 Fisher's Exact Test for Count Data
 
 data:  matrix(10 * 1:9, nc = 3) 
 p-value = 0.3149
 alternative hypothesis: two.sided 
 
 Warning message:
 'hybrid' is ignored for a 2 x 2 table in: fisher.test(matrix(10 * 1:9,
 nc = 3), hybrid = TRUE, workspace = 1e+06) 
 
 whereas (notice the slightly different p value)
 
 .
 .
 .
 

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

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


[R] Forcing the intercept

2006-11-13 Thread Heather Maughan
Dear R-users:

I am doing multiple regressions using the lm function and would like to
force the intercept to be equal to a specific value (such as 4.3).  I was
able to find out how to force it through the origin but this does not work
for other values. 

I am also interested in forcing the regression parameters obtained from one
regression in another regression with a subset of the data.

Are either of these possible in R?  I have been searching the help guide for
hours and have been unsuccessful.

Many thanks,
Heather

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Question on applying vectors to data.frames by row

2006-11-13 Thread RDM
# Newbie alert
# I am wanting to multiply the rows in a dataframe by a vector.
# However, the default behavior appears to be for the vector to be applied
# column wise. For example:

vct - 1:4
df - data.frame(c1 = 5:10, c2= 6:11, c3=7:12, c4=8:13)
multTheTwo - vct * df
multTheTwo

# This results in the vector getting cycled columnwise
#   c1 c2 c3 c4
# 1  5 18  7 24
# 2 12 28 16 36
# 3 21  8 27 10
# 4 32 18 40 22
# 5  9 30 11 36
# 6 20 44 24 52

# But what I actually want is:
#  c1 c2 c3 c4
#1  5 12 21 32which is 5*1, 6*2, 7*3, 8*4
#2  6 14 24 36same pattern applied to the next row
#3  7 16 27 40so on 
#4  8 18 30 44
#5  9 20 33 48
#6 10 22 36 52

# I am currently using a for statement to do this, but that just doesn't
feel
# right. Is there another option that is more R like.

for(i in 1:length(vct)) multTheTwo[,i] - vct[i] * df[,i]
multTheTwo

# Thanks Richard

--

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Profile confidence intervals and LR chi-square test

2006-11-13 Thread Inman, Brant A. M.D.

System: R 2.3.1 on Windows XP machine.

I am building a logistic regression model for a sample of 100 cases in
dataframe d, in which there are 3 binary covariates: x1, x2 and x3.



 summary(d)
 y  x1 x2 x3
 0:54   0:50   0:64   0:78  
 1:46   1:50   1:36   1:22  

 fit - glm(y ~ x1 + x2 + x3, data=d, family=binomial(link=logit))

 summary(fit)

Call:
glm(formula = y ~ x1 + x2 + x3, family = binomial(link = logit), 
data = d)

Deviance Residuals: 
Min   1Q   Median   3Q  Max  
-1.6503  -1.0220  -0.7284   0.9965   1.7069  

Coefficients:
Estimate Std. Error z value Pr(|z|)  
(Intercept)  -0.3772 0.3721  -1.014   0.3107  
x11  -0.8144 0.4422  -1.842   0.0655 .
x21   0.9226 0.4609   2.002   0.0453 *
x31   1.3347 0.5576   2.394   0.0167 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 137.99  on 99  degrees of freedom
Residual deviance: 120.65  on 96  degrees of freedom
AIC: 128.65

Number of Fisher Scoring iterations: 4

 exp(fit$coef)
(Intercept) x11 x21 x31 
  0.6858006   0.4429233   2.5157321   3.7989873 
---

After reading the appropriate sections in MASS4 (7.2 and 8.4 in
particular), I decided to estimate the 95% confidence intervals for the
odds ratios using the profile method implemented in the confint
function. I then used the anova function to perform the deviance
chi-square tests for each covariate.

---
 ci - confint(fit); exp(ci)
Waiting for profiling to be done...
2.5 %97.5 %
(Intercept) 0.3246680  1.413684
x11 0.1834819  1.048154
x21 1.0256096  6.314473
x31 1.3221533 12.129210

 anova(fit, test='Chisq')
Analysis of Deviance Table

Model: binomial, link: logit

Response: y

Terms added sequentially (first to last)


 Df Deviance Resid. Df Resid. Dev P(|Chi|)
NULL99137.989  
x115.85698132.133 0.016
x215.27197126.862 0.022
x316.21296120.650 0.013


My question relates to the interpretation of the significance of
variable x1.  The OR for x1 is 0.443 and its profile confidence interval
is 0.183-1.048.  If a type I error rate of 5% is assumed, this result
would tend to suggest that x1 is NOT a significant predictor of y.
However, the deviance chi-square test has a P-value of 0.016, which
suggests that x1 is indeed a significant predictor of y. How do I
reconcile these two differing messages? I do recognize that the upper
bound of the confidence interval is pretty close to 1, but I am certain
that some journal reviewer will point out the problem as inconsistent.

Brant Inman

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Forcing the intercept

2006-11-13 Thread Leeds, Mark \(IED\)
you can just subtract 4.3 from the independent variable and then do
through zero. That will
Give you a force through 4.3. I don't undersarand the second part of
your statement.



-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Heather Maughan
Sent: Monday, November 13, 2006 6:31 PM
To: r-help@stat.math.ethz.ch
Subject: [R] Forcing the intercept

Dear R-users:

I am doing multiple regressions using the lm function and would like
to force the intercept to be equal to a specific value (such as 4.3).  I
was able to find out how to force it through the origin but this does
not work for other values. 

I am also interested in forcing the regression parameters obtained from
one regression in another regression with a subset of the data.

Are either of these possible in R?  I have been searching the help guide
for hours and have been unsuccessful.

Many thanks,
Heather

__
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
and provide commented, minimal, self-contained, reproducible code.


This is not an offer (or solicitation of an offer) to buy/se...{{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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Forcing the intercept

2006-11-13 Thread Jeffrey Robert Spies
Second part:  you want to look at the residuals after specifying a  
model?  Just use the dependent variables in your subsample as the  
dependent variables in your regression equation and subtract that  
from your outcome variable in  your  subsample.  Might not be the  
answer to the question you're asking though.

Jeff.
-- 
Jeffrey R. Spies
http://www.nd.edu/~jspies/


On Nov 13, 2006, at 7:00 PM, Leeds, Mark (IED) wrote:

 you can just subtract 4.3 from the independent variable and then do
 through zero. That will
 Give you a force through 4.3. I don't undersarand the second part of
 your statement.



 -Original Message-
 From: [EMAIL PROTECTED]
 [mailto:[EMAIL PROTECTED] On Behalf Of Heather Maughan
 Sent: Monday, November 13, 2006 6:31 PM
 To: r-help@stat.math.ethz.ch
 Subject: [R] Forcing the intercept

 Dear R-users:

 I am doing multiple regressions using the lm function and would like
 to force the intercept to be equal to a specific value (such as  
 4.3).  I
 was able to find out how to force it through the origin but this does
 not work for other values.

 I am also interested in forcing the regression parameters obtained  
 from
 one regression in another regression with a subset of the data.

 Are either of these possible in R?  I have been searching the help  
 guide
 for hours and have been unsuccessful.

 Many thanks,
 Heather

 __
 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
 and provide commented, minimal, self-contained, reproducible code.
 

 This is not an offer (or solicitation of an offer) to buy/se... 
 {{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
 and provide commented, minimal, self-contained, reproducible code.


[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Question on applying vectors to data.frames by row

2006-11-13 Thread Gabor Grothendieck
Here are a couple of possibilities:

   as.data.frame(t(t(df) * vct))

   df * rep(vct, each = nrow(df))

On 11/13/06, RDM [EMAIL PROTECTED] wrote:
 # Newbie alert
 # I am wanting to multiply the rows in a dataframe by a vector.
 # However, the default behavior appears to be for the vector to be applied
 # column wise. For example:

 vct - 1:4
 df - data.frame(c1 = 5:10, c2= 6:11, c3=7:12, c4=8:13)
 multTheTwo - vct * df
 multTheTwo

 # This results in the vector getting cycled columnwise
 #   c1 c2 c3 c4
 # 1  5 18  7 24
 # 2 12 28 16 36
 # 3 21  8 27 10
 # 4 32 18 40 22
 # 5  9 30 11 36
 # 6 20 44 24 52

 # But what I actually want is:
 #  c1 c2 c3 c4
 #1  5 12 21 32which is 5*1, 6*2, 7*3, 8*4
 #2  6 14 24 36same pattern applied to the next row
 #3  7 16 27 40so on 
 #4  8 18 30 44
 #5  9 20 33 48
 #6 10 22 36 52

 # I am currently using a for statement to do this, but that just doesn't
 feel
 # right. Is there another option that is more R like.

 for(i in 1:length(vct)) multTheTwo[,i] - vct[i] * df[,i]
 multTheTwo

 # Thanks Richard

 --

 __
 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
 and provide commented, minimal, self-contained, reproducible code.


__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Confidence interval for relative risk

2006-11-13 Thread David Duffy
Michael Dewey [EMAIL PROTECTED] wrote
 Subject: [R] Confidence interval for relative risk

 The concrete problem is that I am refereeing
 a paper where a confidence interval is
 presented for the risk ratio and I do not find
 it credible. I show below my attempts to
 do this in R. The example is slightly changed
 from the authors'.

 I can obtain a confidence interval for
 the odds ratio from fisher.test of
 course

 === fisher.test example ===

  outcome - matrix(c(500, 0, 500, 8), ncol = 2, byrow = TRUE)
  fisher.test(outcome)

  Fisher's Exact Test for Count Data

 data:  outcome
 p-value = 0.00761
 alternative hypothesis: true odds ratio is not equal to 1
 95 percent confidence interval:
   1.694792  Inf
 sample estimates:
 odds ratio
 Inf

 === end example ===


Since

RR   = p1 +1   p1=a/(a+b), p2=c/(c+d)
 --
   p2 1
 -   +   ---
  1-p2OR

you can backsolve for the RR CI.

 x - function(p1, p2, oddsr) p1 + 1/(p2/(1-p2) + 1/oddsr)
 1/x(8/508, 0, 1/1.694792)
[1] 1.650734

Another approach is the Cornfield method -- which is a profile likelihood based 
CI 
using the Gibbs Chi-square (LRTS), but IIRC originally used the Pearson 
chi-square.
I don't think there are R implementations.

David Duffy.

-- 
| David Duffy (MBBS PhD) ,-_|\
| email: [EMAIL PROTECTED]  ph: INT+61+7+3362-0217 fax: -0101  / *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Generating confusion matrix / Kappa stats

2006-11-13 Thread Milton Cezar Ribeiro
Hi there, 
   
  I classified an image and checked out it on the field. Now I have a table 
with three fields like:
   
   
  Field_ID   Field_Class Image_Class
  1 1 1
  2 3 5
  3 4 1
  4 1 1
  5 2 1 
  ...
   
  And now I need gerating a confusion matrix to compute Kappa statistic. First 
of all, how can I generate the confusion matrix from my input table? What 
package is good for compute Kappa statistics?
   
  Regards,
   
  Miltinho
  Brazil
   




-

[[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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Profile confidence intervals and LR chi-square test

2006-11-13 Thread Henric Nilsson
On 2006-11-14 00:41, Inman, Brant A. M.D. skrev:
 System: R 2.3.1 on Windows XP machine.

Time to upgrade!

 
 I am building a logistic regression model for a sample of 100 cases in
 dataframe d, in which there are 3 binary covariates: x1, x2 and x3.

Please provide a reproducible example (as suggested by the posting guide).

 
 
 
 summary(d)
  y  x1 x2 x3
  0:54   0:50   0:64   0:78  
  1:46   1:50   1:36   1:22  
 
 fit - glm(y ~ x1 + x2 + x3, data=d, family=binomial(link=logit))
 
 summary(fit)
 
 Call:
 glm(formula = y ~ x1 + x2 + x3, family = binomial(link = logit), 
 data = d)
 
 Deviance Residuals: 
 Min   1Q   Median   3Q  Max  
 -1.6503  -1.0220  -0.7284   0.9965   1.7069  
 
 Coefficients:
 Estimate Std. Error z value Pr(|z|)  
 (Intercept)  -0.3772 0.3721  -1.014   0.3107  
 x11  -0.8144 0.4422  -1.842   0.0655 .
 x21   0.9226 0.4609   2.002   0.0453 *
 x31   1.3347 0.5576   2.394   0.0167 *
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
 
 (Dispersion parameter for binomial family taken to be 1)
 
 Null deviance: 137.99  on 99  degrees of freedom
 Residual deviance: 120.65  on 96  degrees of freedom
 AIC: 128.65
 
 Number of Fisher Scoring iterations: 4
 
 exp(fit$coef)
 (Intercept) x11 x21 x31 
   0.6858006   0.4429233   2.5157321   3.7989873 
 ---
 
 After reading the appropriate sections in MASS4 (7.2 and 8.4 in
 particular), I decided to estimate the 95% confidence intervals for the
 odds ratios using the profile method implemented in the confint
 function. I then used the anova function to perform the deviance
 chi-square tests for each covariate.
 
 ---
 ci - confint(fit); exp(ci)
 Waiting for profiling to be done...
 2.5 %97.5 %
 (Intercept) 0.3246680  1.413684
 x11 0.1834819  1.048154
 x21 1.0256096  6.314473
 x31 1.3221533 12.129210
 
 anova(fit, test='Chisq')
 Analysis of Deviance Table
 
 Model: binomial, link: logit
 
 Response: y
 
 Terms added sequentially (first to last)
   
Hence, your use of the `anova' function doesn't return tests 
corresponding to the CIs computed above.

 
 
  Df Deviance Resid. Df Resid. Dev P(|Chi|)
 NULL99137.989  
 x115.85698132.133 0.016
 x215.27197126.862 0.022
 x316.21296120.650 0.013
 
 
 My question relates to the interpretation of the significance of
 variable x1.  The OR for x1 is 0.443 and its profile confidence interval
 is 0.183-1.048.  If a type I error rate of 5% is assumed, this result
 would tend to suggest that x1 is NOT a significant predictor of y.

This is also suggested by the Wald test computed by the `summary' function.

 However, the deviance chi-square test has a P-value of 0.016, which
 suggests that x1 is indeed a significant predictor of y. How do I

That p-value corresponds to adding x1 to a model containing only the 
intercept term.

 reconcile these two differing messages? I do recognize that the upper

Generally, in order to compute the LR test for the null hypothesis of 
some subset of the parameters being equal to zero, you need to 
explicitly fit both the restricted and the unrestricted model and 
compare them using the `anova' function.

Also, see FAQ 7.18.


HTH,
Henric



 bound of the confidence interval is pretty close to 1, but I am certain
 that some journal reviewer will point out the problem as inconsistent.
 
 Brant Inman
 
 __
 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
 and provide commented, minimal, self-contained, reproducible code.


__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Correspondence Analysis

2006-11-13 Thread Kris Lockyear
Dear All,

I am in the process of teaching myself R and am getting the hang of it 
slowly, and so apologies for what may be a novice question. (Many thanks to 
those who have helped me so far).

I have been performing normal Correspondence Analysis for a number of 
years using a variety of packages. CA is available in a number of R packages 
(e.g. vegan).

Could anyone advise me which one(s) would provide me with the standard 
Greenacre diagnostic statistics (ie. quality, mass, contribution etc as 
discussed in Correspondence Analysis in Practice, esp. chapter 11) and how I 
get hold of them?

Many thanks in advance, Kris Lockyear.

_
Be the first to hear what's new at MSN - sign up to our free newsletters!

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Building R from source

2006-11-13 Thread YONGWAN CHUN
Hello,


I was trying to build R from source on Windows XP. I installed software which 
are mentioned from the follow web page 
http://www.murdoch-sutherland.com/Rtools/ (Last accessed on Nov. 13th, 2006) . 
Unfortunately, I got error messages whenever I run 'make all recommended' 
without modifying 'MkRules' file. I have removed software and reinstalled them 
several times but I still failed to get it done. The below message is what I 
got.

If anybody gives information, I would appreciate it very much. 

Thanks, 

Yongwan 

===
D:\Rsource\R-2.4.0\src\gnuwin32make all recommended
make[1]: `Rpwd.exe' is up to date.
make[4]: Nothing to be done for `svnonly'.
installing C headers
make[2]: `all' is up to date.
make[2]: `libRblas.dll.a' is up to date.
make[5]: Nothing to be done for `svnonly'.
installing C headers
make --no-print-directory -C ../extra/intl -f Makefile.win
make --no-print-directory -C ../appl  OPTFLAGS='-O3 -Wall -pedantic -std=gnu99' 
FOPTFLAGS='-O3 -Wall' -f Makefile.win
make --no-print-directory -C ../nmath OPTFLAGS='-O3 -Wall -pedantic -std=gnu99' 
-f Makefile.win
make --no-print-directory -C ../main  OPTFLAGS='-O3 -Wall -pedantic -std=gnu99 
-DLEA_MALLOC' FFLAGS='-O3 -Wall' -f Makef
ile.win
make --no-print-directory -C ./graphapp OPTFLAGS='-O3 -Wall -pedantic 
-std=gnu99'
make --no-print-directory -C ./getline OPTFLAGS='-O3 -Wall -pedantic -std=gnu99'
make[4]: `gl.a' is up to date.
make -f Makefile.win chartables.h
make[5]: `chartables.h' is up to date.
make -f Makefile.win makeMakedeps
make -f Makefile.win libpcre.a
make[5]: `libpcre.a' is up to date.
make[4]: Nothing to be done for `all'.
make[4]: Nothing to be done for `all'.
gcc  -shared -s -mwindows -o R.dll R.def console.o dataentry.o dynload.o edit.o 
editor.o embeddedR.o extra.o opt.o pager
.o preferences.o psignal.o rhome.o rui.o run.o shext.o sys-win32.o system.o 
e_pow.o malloc.o ../main/libmain.a ../appl/l
ibappl.a ../nmath/libnmath.a graphapp/ga.a getline/gl.a ../extra/xdr/libxdr.a 
../extra/zlib/libz.a ../extra/pcre/libpcre
.a ../extra/bzip2/libbz2.a ../extra/intl/libintl.a ../extra/trio/libtrio.a 
dllversion.o -L. -lg2c -lRblas -lcomctl32 -lv
ersion
../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x15f9): 
undefined reference to `_pcre_ucp_findprop'
../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1639): 
undefined reference to `_pcre_ucp_findprop'
../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1736): 
undefined reference to `_pcre_ucp_findprop'
../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1770): 
undefined reference to `_pcre_ucp_findprop'
../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x184d): 
undefined reference to `_pcre_ucp_findprop'
../extra/pcre/libpcre.a(pcre_dfa_exec.o):pcre_dfa_exec.c:(.text+0x1e0f): more 
undefined references to `_pcre_ucp_findpro
p' follow
collect2: ld returned 1 exit status
make[3]: *** [R.dll] Error 1
make[2]: *** [../../bin/R.dll] Error 2
make[1]: *** [rbuild] Error 2
make: *** [all] Error 2

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


[R] Problem with file size

2006-11-13 Thread Benilton Carvalho
Hi everyone,

I have 2 environments (2 different R sessions) as described below:

Session 1:

Name of the environment: CrlmmInfo
Objects in the environment:
index1: logical index - length 238304
index2: logical index - length 238304
priors: list of 4 - (matrix 6x6, 2 vectors of length 6, vector of  
length 2) - all num
params: list of 4:
  centers [238304 x 3 x 2]: num
  scales [238304 x 3 x 2]: num
  N [238304 x 3]: num
  f0 [scalar]: num

If I save this environment to a file, I get a file of 23MB. Great.

Session 2:
 Analogous to Session 1, but replace 238304 by 262264.

If I save the environment on Session 2, I get a file of 8.4GB.

I applied object.size on each of the objects in each environment, and  
this is what I got:

For Session 1:
index1: 16204864
index2: 16204864
priors: 3336
params: 74353584

For Session 2:
index1: 1049096
index2: 1049096
priors: 3336
params: 81829104

Is this increase from 23MB to 8.4GB expected to happen?

Benilton

sessionInfo() on both sessions is identical:

  sessionInfo()
R version 2.4.0 (2006-10-03)
x86_64-unknown-linux-gnu

locale:
LC_CTYPE=en_US.iso885915;LC_NUMERIC=C;LC_TIME=en_US.iso885915;LC_COLLATE 
=en_US.iso885915;LC_MONETARY=en_US.iso885915;LC_MESSAGES=en_US.iso885915 
;LC_PAPER=en_US.iso885915;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASU 
REMENT=en_US.iso885915;LC_IDENTIFICATION=C

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

  version
_
platform   x86_64-unknown-linux-gnu
arch   x86_64
os linux-gnu
system x86_64, linux-gnu
status
major  2
minor  4.0
year   2006
month  10
day03
svn rev39566
language   R
version.string R version 2.4.0 (2006-10-03)

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] Installing package rpvm under Windows

2006-11-13 Thread Adrian Dragulescu

Hello,

I'm trying to install the rpvm package under Windows, but I am having
problems.  I have pvm3.4 installed properly.

I've defined the system variables
  PVM_ROOT = C:\PROGRA~1\pvm3.4\
  PVM_ARCH = win32

When I try to install, I get this:
C:\R\PackagesRcmd INSTALL rpvm_1.0.1.tar.gz


-- Making package rpvm 

   **
   WARNING: this package has a configure script
 It probably needs manual configuration
   **

  adding build stamp to DESCRIPTION
  making DLL ...
gcc  -shared -s  -o rpvm.dll rpvm.def rpvm_core.o rpvm_ser.o utils.orpvm_res.o
 -Lc:/PROGRA~1/R/R-24~1.0/bin   -lR
rpvm_core.o:rpvm_core.c:(.text+0x44): undefined reference to `pvm_perror'
rpvm_core.o:rpvm_core.c:(.text+0x49): undefined reference to `pvm_exit'
rpvm_core.o:rpvm_core.c:(.text+0x197): undefined reference to `pvm_mytid'
rpvm_core.o:rpvm_core.c:(.text+0x1c7): undefined reference to `pvm_parent'
rpvm_core.o:rpvm_core.c:(.text+0x207): undefined reference to `pvm_exit'
rpvm_core.o:rpvm_core.c:(.text+0x2a8): undefined reference to `pvm_pstat'
rpvm_core.o:rpvm_core.c:(.text+0x352): undefined reference to `pvm_kill'
rpvm_core.o:rpvm_core.c:(.text+0x3ba): undefined reference to `pvm_tasks'
rpvm_core.o:rpvm_core.c:(.text+0x6bf): undefined reference to `pvm_spawn'
rpvm_core.o:rpvm_core.c:(.text+0x7b7): undefined reference to `pvm_initsend'
rpvm_core.o:rpvm_core.c:(.text+0x7f7): undefined reference to `pvm_mkbuf'
rpvm_core.o:rpvm_core.c:(.text+0x837): undefined reference to `pvm_freebuf'
rpvm_core.o:rpvm_core.c:(.text+0x867): undefined reference to `pvm_getsbuf'
rpvm_core.o:rpvm_core.c:(.text+0x897): undefined reference to `pvm_getrbuf'
rpvm_core.o:rpvm_core.c:(.text+0x8d7): undefined reference to `pvm_setsbuf'
rpvm_core.o:rpvm_core.c:(.text+0x917): undefined reference to `pvm_setrbuf'
rpvm_core.o:rpvm_core.c:(.text+0x97d): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0x9f2): undefined reference to `pvm_pkdouble'
rpvm_core.o:rpvm_core.c:(.text+0xa43): undefined reference to `pvm_pkstr'
rpvm_core.o:rpvm_core.c:(.text+0xa9e): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xad4): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xb2e): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xb5f): undefined reference to `pvm_pkdouble'
rpvm_core.o:rpvm_core.c:(.text+0xbc1): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xbf8): undefined reference to `pvm_pkstr'
rpvm_core.o:rpvm_core.c:(.text+0xcc5): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xcf5): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xdc5): undefined reference to `pvm_pkint'
rpvm_core.o:rpvm_core.c:(.text+0xdf5): undefined reference to `pvm_pkdouble'
rpvm_core.o:rpvm_core.c:(.text+0xe6e): undefined reference to `pvm_send'
rpvm_core.o:rpvm_core.c:(.text+0xee4): undefined reference to `pvm_mcast'
rpvm_core.o:rpvm_core.c:(.text+0xf43): undefined reference to `pvm_recv'
rpvm_core.o:rpvm_core.c:(.text+0xf9e): undefined reference to `pvm_nrecv'
rpvm_core.o:rpvm_core.c:(.text+0xffe): undefined reference to `pvm_probe'
rpvm_core.o:rpvm_core.c:(.text+0x107c): undefined reference to `pvm_trecv'
rpvm_core.o:rpvm_core.c:(.text+0x113d): undefined reference to `pvm_bufinfo'
rpvm_core.o:rpvm_core.c:(.text+0x1217): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x12ac): undefined reference to
`pvm_upkdouble'
rpvm_core.o:rpvm_core.c:(.text+0x130b): undefined reference to `pvm_upkstr'
rpvm_core.o:rpvm_core.c:(.text+0x1374): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x13b4): undefined reference to `pvm_upkstr'
rpvm_core.o:rpvm_core.c:(.text+0x1421): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x146f): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x14c1): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x150f): undefined reference to
`pvm_upkdouble'
rpvm_core.o:rpvm_core.c:(.text+0x1561): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x15be): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x1611): undefined reference to `pvm_upkint'
rpvm_core.o:rpvm_core.c:(.text+0x166e): undefined reference to
`pvm_upkdouble'
rpvm_core.o:rpvm_core.c:(.text+0x16c0): undefined reference to `pvm_config'
rpvm_core.o:rpvm_core.c:(.text+0x1894): undefined reference to
`pvm_start_pvmd'
rpvm_core.o:rpvm_core.c:(.text+0x191c): undefined reference to
`pvm_addhosts'
rpvm_core.o:rpvm_core.c:(.text+0x19c7): undefined reference to
`pvm_delhosts'
rpvm_core.o:rpvm_core.c:(.text+0x1a27): undefined reference to `pvm_halt'
rpvm_core.o:rpvm_core.c:(.text+0x1acc): undefined reference to `pvm_mstat'
rpvm_core.o:rpvm_core.c:(.text+0x1b53): undefined reference to
`pvm_joingroup'
rpvm_core.o:rpvm_core.c:(.text+0x1ba3): undefined reference to 

Re: [R] Using lrm

2006-11-13 Thread Andrew Robinson
Hello Nitin,

if you examine the help information for lrm carefully, at the bottom
you will find numerous examples that you can follow.

?lrm

By the way, asking a queestion like this it's best to clarify if you
mean the lrm from the Design package or some other one.  I assume it's
from Design, as the syntax you quote is a good match.

It's also best to study the help files carefully.

Andrew

On Mon, Nov 13, 2006 at 11:01:59PM -0600, nitin jindal wrote:
 Hi,
 
 I have to build a logistic regression model on a data set that I have. I
 have three input variables (x1, x2, x3) and one output variable (y).
 
 The syntax of lrm function looks like this
 
  lrm(formula, data, subset, na.action=na.delete, method=lrm.fit,
  model=FALSE, x=FALSE, y=FALSE, linear.predictors=TRUE, se.fit=FALSE,
 
  penalty=0, penalty.matrix, tol=1e-7,
  strata.penalty=0, var.penalty=c('simple','sandwich'),
  weights, normwt, ...)
 
 Any logistic regression model would take y and x1,x2,x3 as parameters and
 output the model (probabilities). So, I dont know where to fit in these
 values in this function.
 
 Any help is appreciated. I am chasing a deadline in my project.
 
 Nitin
 
   [[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
 and provide commented, minimal, self-contained, reproducible code.

-- 
Andrew Robinson  
Department of Mathematics and StatisticsTel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

__
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
and provide commented, minimal, self-contained, reproducible code.


[R] 回覆: Re: Need help in waveslim package: imodwt and universal.thresh.modwt

2006-11-13 Thread Airon Yiu
Hi  Rogerio:
   
  I am using Waveslim 1.5 on R 2.3.0, running on Chinese WinXP (SP2).  
   
  The data, attached in the CSV file, is a stock data (0001 from Hong Kong) 
downloaded from Yahoo!Finance.  The time series data are the Adj. Close prices 
(last column) from 4-Jan-00 to 30-Nov-2005.
   
  The waveslim mra rountines are working fine on these data.
   
  Thks

rdporto1 [EMAIL PROTECTED] 說:
  The first two commands worked fine to me. I used
xdata = rnorm(128).

To let us help you more, specify your R version,
some data, OS etc, as usually asked.

Rogerio Porto.

-- Cabe蓷lho original ---

De: [EMAIL PROTECTED]
Para: r-help@stat.math.ethz.ch
C鏕ia:
Data: Sun, 12 Nov 2006 16:02:45 +0800 (CST)
Assunto: [R] Need help in waveslim package: imodwt and universal.thresh.modwt

 Hi:
 I have encountered problems with imodwt and universal.thresh.modwt and cannot 
 find any reference in R Search. Hope someone can give me some ideas:

 Starting with
 modwt.la8 - modwt(xdata, la8, n.level=6) -- this seems to work fine

 (1) ydata - imodwt(modwt.la8)
 will always give ydata as numeric(0) (no values) instead of being a time 
 series data with the same lenght as xdata.

 (2) thred.la8 - universal.thresh.modwt(modwt.la8, max.level=4, hard=F) will 
 give me error at:
 abs(wc.fine) - cannot contain non-numeric field.

 Any help is much appreciated.

 Regards







Content-Type: text/html; charset=big5
Content-Transfer-Encoding: 8bit

DIVHinbsp; Rogerio:/DIV  DIVnbsp;/DIV  DIVI am using Waveslim 1.5 
on R 2.3.0, running on Chinese WinXP (SP2).nbsp; /DIV  DIVnbsp;/DIV  
DIVThe data, attached in the CSV file, is a stock data (0001 from Hong Kong) 
downloaded from Yahoo!Finance.nbsp; The time series data are thenbsp;Adj. 
Close prices (last column) fromnbsp;4-Jan-00 to 30-Nov-2005./DIV  
DIVnbsp;/DIV  DIVThe waveslimnbsp;mra rountines are working fine on 
these data./DIV  DIVnbsp;/DIV  DIVThksBRBRBIrdporto1 
lt;[EMAIL PROTECTED]gt;/I/B 說:/DIV  BLOCKQUOTE class=replbq 
style=PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #1010ff 2px solidThe 
first two commands worked fine to me. I usedBRxdata = rnorm(128).BRBRTo 
let us help you more, specify your R version,BRsome data, OS etc, as usually 
asked.BRBRRogerio Porto.BRBR-- Cabe蓷lho original 
---BRBRDe: [EMAIL PROTECTED]BRPara:
 r-help@stat.math.ethz.chBRC鏕ia:BRData: Sun, 12 Nov 2006 16:02:45 +0800 
(CST)BRAssunto: [R] Need help in waveslim package: imodwt and 
universal.thresh.modwtBRBRgt; Hi:BRgt; I have encountered problems with 
imodwt and universal.thresh.modwt and cannot find any reference in R Search. 
Hope someone can give me some ideas:BRgt;BRgt; Starting withBRgt; 
modwt.la8 lt;- modwt(xdata, la8, n.level=6) lt;-- this seems to work 
fineBRgt;BRgt; (1) ydata lt;- imodwt(modwt.la8)BRgt; will always give 
ydata as numeric(0) (no values) instead of being a time series data with the 
same lenght as xdata.BRgt;BRgt; (2) thred.la8 lt;- 
universal.thresh.modwt(modwt.la8, max.level=4, hard=F) will give me error 
at:BRgt; abs(wc.fine) - cannot contain non-numeric field.BRgt;BRgt; 
Any help is much appreciated.BRgt;BRgt; 
RegardsBRgt;BRgt;BRBR/BLOCKQUOTEBRp#32;___br
 YM - 離線訊息br

Content-Type: application/octet-stream; name=CKH-0001-ALL-data-only.csv
Content-Transfer-Encoding: base64
Content-Description: 3176266068-CKH-0001-ALL-data-only.csv
Content-Disposition: attachment; filename=CKH-0001-ALL-data-only.csv

MDg2MTEsNzkuNjUNCg==

__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Problem with file size

2006-11-13 Thread Prof Brian Ripley
On Mon, 13 Nov 2006, Benilton Carvalho wrote:

 Hi everyone,

 I have 2 environments (2 different R sessions) as described below:

 Session 1:

 Name of the environment: CrlmmInfo
 Objects in the environment:
index1: logical index - length 238304
index2: logical index - length 238304
priors: list of 4 - (matrix 6x6, 2 vectors of length 6, vector of
 length 2) - all num
params: list of 4:
  centers [238304 x 3 x 2]: num
  scales [238304 x 3 x 2]: num
  N [238304 x 3]: num
  f0 [scalar]: num

 If I save this environment to a file, I get a file of 23MB. Great.

 Session 2:
 Analogous to Session 1, but replace 238304 by 262264.

 If I save the environment on Session 2, I get a file of 8.4GB.

 I applied object.size on each of the objects in each environment, and
 this is what I got:

 For Session 1:
 index1: 16204864
 index2: 16204864
 priors: 3336
 params: 74353584

 For Session 2:
 index1: 1049096
 index2: 1049096
 priors: 3336
 params: 81829104

 Is this increase from 23MB to 8.4GB expected to happen?

We don't have enough information to expect anything.  Saving environments 
is a tricky business: read the description in the R Internals manual if 
this is news to you.

But note that the index[12] objects are much smaller in Session 2, and 
those in Session 1 are about 16x larger than needed for the description 
given.  So my guess is that they have attributes which are promises in 
Session 2, and save() there is pulling in another environment to enable 
those attributes to be evaluated later.

-- 
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Embedded carriage returns in text document

2006-11-13 Thread Jeffrey Robert Spies
For OS X, http://dos2unix.darwinports.com/

--  
Jeffrey R. Spies
http://www.nd.edu/~jspies/


On Nov 13, 2006, at 11:06 PM, [EMAIL PROTECTED]  
[EMAIL PROTECTED] wrote:

 Dennis,

 You can get rid of the '^M' by reformating the
 file from DOS to UNIX. Withing a UNIX system
 send the command
 dos2unix filename1 filename2.
 Filename2 won't have the '^M'.

 Hope it helps,

 Augusto


 
 Augusto Sanabria. MSc, PhD.
 Mathematical Modeller
 Risk Research Group
 Geospatial  Earth Monitoring Division
 Geoscience Australia (www.ga.gov.au)
 Cnr. Jerrabomberra Av.  Hindmarsh Dr.
 Symonston ACT 2601
 Ph. (02) 6249-9155






 -Original Message-
 From: [EMAIL PROTECTED]
 [mailto:[EMAIL PROTECTED] On Behalf Of Dennis Fisher
 Sent: Tuesday, 14 November 2006 9:06
 To: r-help@stat.math.ethz.ch
 Subject: [R] Embedded carriage returns in text document


 Colleagues,

 I am using R 2.4.0 on both a Mac (10.4.8) and Linux (RedHat 9).  To
 read data from an Excel spreadsheet, I do save as  in Excel, then
 select the Text (tab-delimited) format.  The resulting file uses a
 tab separator and I can usually read the file using read.delim.

 Sometimes, the header row contains embedded carriage returns.  When I
 view the file, these carriage returns appear as  ^M.

 Now the problem:
 When I read.delim these files, they do not read correctly.  Sometimes
 I get error messages; sometimes only the first line is read.
 Interestingly, invoking the option skip=1 (or a larger N) does not
 appear to bypass the problem.

 I can solve the problem by manually deleting these carriage returns
 either in the original Excel file or the .txt version.   However,
 this is not an ideal solution.

 Does anyone have a work-around within R?

 Dennis


 Dennis Fisher MD
 P  (The P Less Than Company)
 Phone: 1-866-PLessThan (1-866-753-7784)
 Fax: 1-415-564-2220
 www.PLessThan.com 


   [[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
 and provide commented, minimal, self-contained, reproducible code.

 __
 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
 and provide commented, minimal, self-contained, reproducible code.


[[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
and provide commented, minimal, self-contained, reproducible code.


[R] Variance of a complex estimator using survey package ...

2006-11-13 Thread justin bem
Hi, 

I want to compute the variance of two complex statistics. The first statistic 
is the a ratio 
   R1=Q(a1)/Q(a2) 
where Q denote de quantile at a1 and a2. The second is also a ratio but not a 
classic one this ratio is 
   R2=sum(x_{i}|x_{i}Q(a1))/sum(x_{i} }|Q(a2))

Linearisation for those statistic is very complex. I want to use bootstrap and 
JKn How can I procede since I dont have a function like svrepanyfunction where 
you just specify the function to compute the statistic.

Sincerly. 


Justin BEM
Elève Ingénieur Statisticien Economiste
BP 294 Yaoundé.
Tél (00237)9597295.






___ 
Découvrez une nouvelle façon d'obtenir des réponses à toutes vos questions ! 
Profitez des connaissances, des opinions et des expériences des internautes sur 
Yahoo! Questions/Réponses 

[[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
and provide commented, minimal, self-contained, reproducible code.


[R] reconstruction of data after PCA and interpolation

2006-11-13 Thread Poizot Emmanuel

Dear all,

I have a matrix of size N x M. I purchase a PCA analysis through prcomp 
function. Then I keep the H eigenvectors which explain 90% of the total 
variance and interpolate each vectors of the matrix H x M, to obtain a 
new matrix of size H x K (K  M). My question is : from this last matrix 
(H x K), how to (re)construct a matrix of size N x K which have the same 
units of the inital matrix (N x M).

Regards

--
Cordialement


Emmanuel Poizot
Cnam/Intechmer
B.P. 324
50103 Cherbourg Cedex

Phone (Direct) : (00 33)(0)233887342
Fax : (00 33)(0)233887339


__
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
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Installing package rpvm under Windows

2006-11-13 Thread Prof Brian Ripley
Plese heed the warning enclosed in ***.  You haven't set up the 
pvm libraries in the link.

If you don't know how to do the manual configuration, talk to the package 
maintainer (as the posting guide asked you to do before posting).

On Mon, 13 Nov 2006, Adrian Dragulescu wrote:

 I'm trying to install the rpvm package under Windows, but I am having
 problems.  I have pvm3.4 installed properly.

 I've defined the system variables
  PVM_ROOT = C:\PROGRA~1\pvm3.4\

Using \ is unlikely to work with the tools that we use, but I see no sign 
that this is being used.

  PVM_ARCH = win32

 When I try to install, I get this:
 C:\R\PackagesRcmd INSTALL rpvm_1.0.1.tar.gz


 -- Making package rpvm 

   **
   WARNING: this package has a configure script
 It probably needs manual configuration
   **

  adding build stamp to DESCRIPTION
  making DLL ...
 gcc  -shared -s  -o rpvm.dll rpvm.def rpvm_core.o rpvm_ser.o utils.orpvm_res.o
 -Lc:/PROGRA~1/R/R-24~1.0/bin   -lR
 rpvm_core.o:rpvm_core.c:(.text+0x44): undefined reference to `pvm_perror'

[...]

-- 
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
and provide commented, minimal, self-contained, reproducible code.