Re: [R] Trellis.par.set/ family/ global change font?

2006-06-30 Thread Renaud Lancelot
A radical solution is to edit the file .../etc/Rdevga and to change
the definition of the first set of fonts, e.g.

# TT Arial : plain
# TT Arial : bold
# TT Arial : italic
# TT Arial : bolditalic

TT Times New Roman : plain
TT Times New Roman : bold
TT Times New Roman : italic
TT Times New Roman : bolditalic

Best,

Renaud

2006/6/30, Deepayan Sarkar [EMAIL PROTECTED]:
 On 6/29/06, McClatchie, Sam (PIRSA-SARDI)
 [EMAIL PROTECTED] wrote:
  Background:
  OS: Linux Ubuntu Dapper 6.06
  release: R 2.3.1
  editor: GNU Emacs 21.4.1
  front-end: ESS 5.2.3
  -
 
  Colleagues
 
  I have a rather complicated trellis plot that a journal editor has 
  requested I edit and change all the fonts to times.
 
  I'd like to change all fonts globally for the plot, as in 
  par(family=serif) for non-trellis plots. Various experiments with 
  trellis.par.set after reading the help page have not solved the problem for 
  me. Even doing the local change
  trellis.par.set(par.xlab.text=list(cex=1.5, family=serif)) does not 
  change the font to times for xlab (I mean my syntax is wrong, not that 
  there is a bug).
 
  So I'm obviously misreading the help page or just missing the meaning. Any 
  suggestions?

 Lattice uses 'fontfamily' rather than 'family' (borrowed from grid, I
 suppose). I don't think there's a way to set the family globally. You
 might try

 trellis.par.set(grid.pars = list(fontfamily = serif))

 but I'm not sure if that will work.

 -Deepayan

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



-- 
Renaud LANCELOT
Département Elevage et Médecine Vétérinaire (EMVT) du CIRAD
Directeur adjoint chargé des affaires scientifiques

CIRAD, Animal Production and Veterinary Medicine Department
Deputy director for scientific affairs

Campus international de Baillarguet
TA 30 / B (Bât. B, Bur. 214)
34398 Montpellier Cedex 5 - France
Tél   +33 (0)4 67 59 37 17
Secr. +33 (0)4 67 59 39 04
Fax   +33 (0)4 67 59 37 95

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


Re: [R] Trellis.par.set/ family/ global change font?

2006-06-30 Thread Prof Brian Ripley
On Fri, 30 Jun 2006, Renaud Lancelot wrote:

 A radical solution is to edit the file .../etc/Rdevga and to change
 the definition of the first set of fonts, e.g.

 # TT Arial : plain
 # TT Arial : bold
 # TT Arial : italic
 # TT Arial : bolditalic

 TT Times New Roman : plain
 TT Times New Roman : bold
 TT Times New Roman : italic
 TT Times New Roman : bolditalic

He is on Linux, so that will not work.  We have not been told the graphics 
device, but both postscript() and pdf() allow the device to be set when 
the device is opened, and that works for me:

pdf(family=serif)
example(xyplot)
dev.off()



 Best,

 Renaud

 2006/6/30, Deepayan Sarkar [EMAIL PROTECTED]:
 On 6/29/06, McClatchie, Sam (PIRSA-SARDI)
 [EMAIL PROTECTED] wrote:
 Background:
 OS: Linux Ubuntu Dapper 6.06
 release: R 2.3.1
 editor: GNU Emacs 21.4.1
 front-end: ESS 5.2.3
 -

 Colleagues

 I have a rather complicated trellis plot that a journal editor has 
 requested I edit and change all the fonts to times.

 I'd like to change all fonts globally for the plot, as in 
 par(family=serif) for non-trellis plots. Various experiments with 
 trellis.par.set after reading the help page have not solved the problem for 
 me. Even doing the local change
 trellis.par.set(par.xlab.text=list(cex=1.5, family=serif)) does not 
 change the font to times for xlab (I mean my syntax is wrong, not that 
 there is a bug).

 So I'm obviously misreading the help page or just missing the meaning. Any 
 suggestions?

 Lattice uses 'fontfamily' rather than 'family' (borrowed from grid, I
 suppose). I don't think there's a way to set the family globally. You
 might try

 trellis.par.set(grid.pars = list(fontfamily = serif))

 but I'm not sure if that will work.

 -Deepayan

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





-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

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


Re: [R] cat and positioning of the output

2006-06-30 Thread J?rn Schulz

Thanks at all for your carefully help. I decided to used both version of your
proposals and to combine with a request of the operating system.

All the best - J?rn Schulz.
-- 
View this message in context: 
http://www.nabble.com/cat-and-positioning-of-the-output-tf1869521.html#a5115737
Sent from the R help forum at Nabble.com.

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


Re: [R] Empirical CDF

2006-06-30 Thread Peter Dalgaard
[EMAIL PROTECTED] writes:

 Good day everyone,
 
 I want to assess the error when fitting a Gram-Charlier
 CDF to some data 'ws', that is, I want to calculate:
 
  Err = |ecdf(ws) - GCh_ser(ws)|
 
 The problem is, I cannot get the F(x) values from the
 ecdf.
 
 'Summary(ecdf())' returns some of the x-axis values,
 but how do you get the F(x) values?

By applying the function returned by ecdf?

Fn - ecdf(some.data)
Fn(some.other.data)

should work. 
 
 Thank you for any help you can provide.
 
 Regards,
 
 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 2609
 Ph. (02) 6249-9155
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
 

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

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


Re: [R] Cointegration Test in R

2006-06-30 Thread Pfaff, Bernhard Dr.
Hello Dennis,

have you considered the function bh6lrtest() in the package urca? 
To my knowledge, there is no other package available that offers VECM 
functionalities. 

Best,
Bernhard

ps: As you migth be interested in VAR and SVAR too: I am currently working on 
such a package which should be submitted to CRAN after summer. 

Dr. Bernhard Pfaff
Global Structured Products Group
(Europe)

Invesco Asset Management Deutschland GmbH
Bleichstrasse 60-62
D-60313 Frankfurt am Main

Tel: +49(0)69 29807 230
Fax: +49(0)69 29807 178
Email: [EMAIL PROTECTED]  

-Ursprüngliche Nachricht-
Von: [EMAIL PROTECTED] 
[mailto:[EMAIL PROTECTED] Im Auftrag von 
Dennis Heidemann
Gesendet: Donnerstag, 29. Juni 2006 21:20
An: R-help@stat.math.ethz.ch
Betreff: [R] Cointegration Test in R

Hello!

I'm using the blrtest() function in the urca package
to test cointegration relationships.
Unfortunately, the hypothesis (restrictions on beta) 
specifies the same restriction on all cointegration vectors.
Is there any possibility to specify different restrictions on
the cointegration vectors?
Are there any other packages in R using cointegration tests?

Thanks and best regards.
Dennis Heidemann
   [[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

*
Confidentiality Note: The information contained in this mess...{{dropped}}

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


Re: [R] R server

2006-06-30 Thread john seers \(IFR\)

There is the Rserve package which you can look at here:

http://stats.math.uni-augsburg.de/Rserve/down.shtml


JS


 
---

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


Re: [R] smoothing splines and degrees of freedom

2006-06-30 Thread i.m.s.white
Steven,

I cannot vouch for the behaviour of the function smooth.spline(), but
the theoretical answer to your question is yes. If g = Sy is the transformation
from data vector y to spline vector g, the equivalent degrees of freedom are
usually defined as EDF = trace(S), where S is the n x n smoothing matrix:

EDF = sum_i(1/(1+theta*lambda_i)),

where lambda_1 to lambda_n are the eigenvalues of S. Two of these are zero, so

EDF = 2 + sum(1/(1+theta*lambda_i))

the sum now over i=3 to n. Here theta is the smoothing parameter. Setting
theta = 0 (no smoothing) gives EDF=n and produces the interpolating spline.
Setting theta = infty gives EDF=2 and a straight line fit. See either

Green and Silverman, Nonparametric regression and generalized linear models,
(p37), or
Hastie and Tibshirani, Generalized additive models, p52.

On Sat, Jun 24, 2006 at 11:35:16AM -0400, Steven Shechter wrote:
 Hi,
 If I set df=2 in my smooth.spline function, is that equivalent to running
 a linear regression through my data?  It appears that df=# of data points
 gives the interpolating spline and that df = 2 gives the linear
 regression, but I just want to confirm this.
 
 Thank you,
 Steven
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 

*I.White   *
*University of Edinburgh   *
*Ashworth Laboratories, West Mains Road*
*Edinburgh EH9 3JT *
*Fax: 0131 650 6564   Tel: 0131 650 5490   *
*E-mail: [EMAIL PROTECTED]  *

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


Re: [R] R server

2006-06-30 Thread roger bos
Mike, may also want to check out Rpad.  Its mainly used to delivery R
programs to other users via a web browser, so the other users don't need R
installed on their machine.  If you make a nice HTML gui, the other users
don't even have to know R.  Each user gets his own environment, so one
person's variables won't interfer with anothers.

HTH,
Roger


On 6/30/06, john seers (IFR) [EMAIL PROTECTED] wrote:


 There is the Rserve package which you can look at here:

 http://stats.math.uni-augsburg.de/Rserve/down.shtml


 JS



 ---

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


[[alternative HTML version deleted]]

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


[R] data extraction

2006-06-30 Thread yohannes alazar
Dear mailing list I have a data that have 20,000 rows and 20 columns. Io
wonted to extract the 10th row only. Example the 10th, 20th, 30th 40th…..2
th. can you please help me how do I do that.Than kyou.
Example is below.
Inpute:
AG GG GG AG
CC CC CC CC
CT CC CT CT
GG GG GG GG
CC CC CC CC
GG GG GG GG
CC CC CC CC
GG CG CG GG
GG GG GG GG
*CC CC CC CC*
AA AG AG AA
AA AA AA AA
GG AG AG GG
GG AG AG GG
GG GG GG GG
TT TT TT TT
AA AG AG AA
CT CC CT TT
AG AA AG GG
*AA AA AG AG*
AA AA AA AA
CC CC CC CG
GG GG GG AG
CT TT CT CT
AT AA AT AT
GG GG GG AG
CG CC CG GG
GG GG GG AG
CC CC CC CT
GG GG GG GG
*GG GG AG AG
*CC CC CC CT
TT TT TT CT
AG GG AG AG
GG GG GG GG
AG AG AA AA
AG GG AG AA
CT TT CT CT
CT CT CC CT
*AC CC AC AC*

output
*CC CC CC CC
AA AA AG AG
GG GG GG GG
AC CC AC AC*

thankyou a lot inadvance!

[[alternative HTML version deleted]]

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

Re: [R] data extraction

2006-06-30 Thread Andris Jankevics
This should work:

data[seq(1,nrow(data),10),]



Andris

On Piektdiena, 30. Jūnijs 2006 14:45, yohannes alazar wrote:
 Dear mailing list I have a data that have 20,000 rows and 20 columns. Io
 wonted to extract the 10th row only. Example the 10th, 20th, 30th
 40th�..2 th. can you please help me how do I do that.Than kyou.
 Example is below.
 Inpute:
 AG GG GG AG
 CC CC CC CC
 CT CC CT CT
 GG GG GG GG
 CC CC CC CC
 GG GG GG GG
 CC CC CC CC
 GG CG CG GG
 GG GG GG GG
 *CC CC CC CC*
 AA AG AG AA
 AA AA AA AA
 GG AG AG GG
 GG AG AG GG
 GG GG GG GG
 TT TT TT TT
 AA AG AG AA
 CT CC CT TT
 AG AA AG GG
 *AA AA AG AG*
 AA AA AA AA
 CC CC CC CG
 GG GG GG AG
 CT TT CT CT
 AT AA AT AT
 GG GG GG AG
 CG CC CG GG
 GG GG GG AG
 CC CC CC CT
 GG GG GG GG
 *GG GG AG AG
 *CC CC CC CT
 TT TT TT CT
 AG GG AG AG
 GG GG GG GG
 AG AG AA AA
 AG GG AG AA
 CT TT CT CT
 CT CT CC CT
 *AC CC AC AC*

 output
 *CC CC CC CC
 AA AA AG AG
 GG GG GG GG
 AC CC AC AC*

 thankyou a lot inadvance!

   [[alternative HTML version deleted]]

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

Re: [R] Extracting R plots from MS Word

2006-06-30 Thread maj
Just a note to say what I did. I think that the results were OK but I have
yet to hear from the journal.

1. I saved the Word document under another name.
2. I deleted all the contents of the document except the target graphic.
3. I printed to file yielding a .prn file.
4. I changed the extension to .ps.
5. I used Gsview PS to EPS.

Murray Jorgensen

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


Re: [R] data extraction

2006-06-30 Thread Chuck Cleland
mydf[seq(10,2,10),]

yohannes alazar wrote:
 Dear mailing list I have a data that have 20,000 rows and 20 columns. Io
 wonted to extract the 10th row only. Example the 10th, 20th, 30th 40th…..2
 th. can you please help me how do I do that.Than kyou.
 Example is below.
 Inpute:
 AG GG GG AG
 CC CC CC CC
 CT CC CT CT
 GG GG GG GG
 CC CC CC CC
 GG GG GG GG
 CC CC CC CC
 GG CG CG GG
 GG GG GG GG
 *CC CC CC CC*
 AA AG AG AA
 AA AA AA AA
 GG AG AG GG
 GG AG AG GG
 GG GG GG GG
 TT TT TT TT
 AA AG AG AA
 CT CC CT TT
 AG AA AG GG
 *AA AA AG AG*
 AA AA AA AA
 CC CC CC CG
 GG GG GG AG
 CT TT CT CT
 AT AA AT AT
 GG GG GG AG
 CG CC CG GG
 GG GG GG AG
 CC CC CC CT
 GG GG GG GG
 *GG GG AG AG
 *CC CC CC CT
 TT TT TT CT
 AG GG AG AG
 GG GG GG GG
 AG AG AA AA
 AG GG AG AA
 CT TT CT CT
 CT CT CC CT
 *AC CC AC AC*
 
 output
 *CC CC CC CC
 AA AA AG AG
 GG GG GG GG
 AC CC AC AC*
 
 thankyou a lot inadvance!
 
   [[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

-- 
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894

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


Re: [R] data extraction

2006-06-30 Thread Dimitris Rizopoulos
you need something like:

new.data - data[seq(10, nrow(data), 10), ]


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


- Original Message - 
From: yohannes alazar [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Sent: Friday, June 30, 2006 1:45 PM
Subject: [R] data extraction


Dear mailing list I have a data that have 20,000 rows and 20 columns. 
Io
wonted to extract the 10th row only. Example the 10th, 20th, 30th 
40th...2
th. can you please help me how do I do that.Than kyou.
Example is below.
Inpute:
AG GG GG AG
CC CC CC CC
CT CC CT CT
GG GG GG GG
CC CC CC CC
GG GG GG GG
CC CC CC CC
GG CG CG GG
GG GG GG GG
*CC CC CC CC*
AA AG AG AA
AA AA AA AA
GG AG AG GG
GG AG AG GG
GG GG GG GG
TT TT TT TT
AA AG AG AA
CT CC CT TT
AG AA AG GG
*AA AA AG AG*
AA AA AA AA
CC CC CC CG
GG GG GG AG
CT TT CT CT
AT AA AT AT
GG GG GG AG
CG CC CG GG
GG GG GG AG
CC CC CC CT
GG GG GG GG
*GG GG AG AG
*CC CC CC CT
TT TT TT CT
AG GG AG AG
GG GG GG GG
AG AG AA AA
AG GG AG AA
CT TT CT CT
CT CT CC CT
*AC CC AC AC*

output
*CC CC CC CC
AA AA AG AG
GG GG GG GG
AC CC AC AC*

thankyou a lot inadvance!

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


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


[R] Biobass, SAGx, and Jonckheere-Terpstra test

2006-06-30 Thread Chris Andrews
Horace,
Last year I found another Jonckheere-Terpstra function at 
http://www.biostat.wustl.edu/archives/html/s-news/2000-10/msg00126.html
Using that and the version in SAGx and wanting a few other options, I 
created the version below.

Hope it is useful,
Chris
-- 
Christopher Andrews, PhD
SUNY Buffalo, Department of Biostatistics



#
#
#
j.test - function(x, y,
   alternative = c(two.sided, decreasing, increasing),
   asymp=TRUE, correct=FALSE, perm=0, na.action=c(omit, fail),
   permgraph=FALSE, permreps=FALSE, ...) {

#Jonckheere-Terpstra test
#x = response
#y = group membership
#alternative hypothesis
#asymptotic = asymptotic formula for variance or don't bother
#correct = continuity correction
#perm = number of repetitions for a permutation test
#na.action = what to do with missing data
#permgraph = draw a histogram of the permutations?
#permreps = return the permutations?
#... for graph

   alternative - match.arg(alternative)
   complete - complete.cases(x,y)
   nn - sum(complete)
   cat(nn, complete cases.\n)
   if (!all(complete)  (na.action==fail)) {
 cat(option na.action is 'fail' and,
   sum(!complete), cases are not complete.\n)
 return(NULL)
   }
   response - x[complete]
   groupvec - y[complete, drop=TRUE]

   if(!is.ordered(groupvec)) {
 groupvec - as.ordered(groupvec)
 cat(y was not an ordered factor.  Redefined to be one.\n)
   }

   lev - levels(groupvec)
   nlev - length(lev)
   if(nlev = 2) {
 if (nlev  2) {
   cat(Not enough groups (k =, nlev, ) for Jonckheere.\n)
   return(NULL)
 } else {
   cat(Two groups.  Could use wilcox.test.\n)
 }
   }

   computestat - function(response, groupvec, n.groups) {
 pairwise - function(x, y) {
   length(x)*(length(y) + (1+length(x))/2) - 
sum(rank(c(x,y))[seq(along=x)])
 }
 H-0
 for (i in seq(n.groups-1))
   for (j in seq(i+1, n.groups))
 H - H+pairwise(response[groupvec == lev[i]], response[groupvec 
== lev[j]])
 return(H)
   }

   retval - list(statistic=computestat(response, groupvec, nlev),
 alternative = paste(alternative,
   paste(levels(groupvec),collapse=switch(alternative,
 two.sided = , , decreasing=   , increasing =   )), 
sep=: ),
 method = Jonckheere,
 data.name = paste(deparse(substitute(x)), by, 
deparse(substitute(y

   if (asymp) {
 ns - tabulate(as.numeric(groupvec))
 retval$EH - (nn * nn - sum(ns * ns))/4
 retval$VH - if (!any(duplicated(response))) {
   (nn * nn * (2 * nn + 3) - sum(ns * ns * (2 * ns + 3)))/72
 } else {
   ds - as.vector(table(response))
   ((nn*(nn-1)*(2*nn+5) - sum(ns*(ns-1)*(2*ns+5)) - 
sum(ds*(ds-1)*(2*ds+5)))/72 +
   sum(ns*(ns-1)*(ns-2))*sum(ds*(ds-1)*(ds-2))/(36*nn*(nn-1)*(nn-2)) +
   sum(ns*(ns-1))*sum(ds*(ds-1))/(8*nn*(nn - 1)))
 }
 pp - if (!correct) {
   pnorm(retval$statistic, retval$EH, sqrt(retval$VH))
 } else if (retval$statistic = retval$EH + 1) {
   pnorm(retval$statistic - 1, retval$EH, sqrt(retval$VH))
 } else if (retval$statistic = retval$EH - 1) {
   pnorm(retval$statistic + 1, retval$EH, sqrt(retval$VH))
 } else {
   .5
 }
 retval$p.value - switch(alternative,
   two.sided = 2*min(pp,1-pp), decreasing = pp, increasing = 1-pp)
   }
   if (perm0) {
 reps - numeric(perm)
 for (i in seq(perm)) {
   reps[i] - computestat(response, sample(groupvec), nlev)
 }
 retval$EH.perm - mean(reps)
 retval$VH.perm - var(reps)
 ppp - if (!correct) {
   pnorm(retval$statistic, retval$EH.perm, sqrt(retval$VH.perm))
 } else if (retval$statistic = retval$EH.perm + 1) {
   pnorm(retval$statistic - 1, retval$EH.perm, sqrt(retval$VH.perm))
 } else if (retval$statistic = retval$EH.perm - 1) {
   pnorm(retval$statistic + 1, retval$EH.perm, sqrt(retval$VH.perm))
 } else {
   .5
 }
 retval$p.value.perm - switch(alternative,
   two.sided = 2*min(ppp,1-ppp), decreasing = ppp, increasing = 1-ppp)
 rr - rank(c(retval$statistic, reps))[1]/(perm+2)
 retval$p.value.rank - switch(alternative,
   two.sided = 2*min(rr,1-rr), decreasing = rr, increasing = 1-rr)
 if (permreps)
   retval$reps - reps
 if (permgraph) {
   hist(reps, xlab=Jonckheere Statistic, prob=TRUE,
 main=Histogram of Permutations, sub=paste(perm, permutations))
   abline(v=retval$statistic, col=red)
   points((-3:3)*sqrt(retval$VH.perm)+retval$EH.perm, rep(0,7), 
col=blue,
 pch=as.character(c(3:1,0:3)))
 }
   }
   class(retval) - htest
   names(retval$statistic) - J
   return(retval)

#returns list (of class htest) with the following components
#Always:
#statistic = value of J
#alternative = same as input
#method = Jonckheere
#data.name = source of data based on call
#Most of the time (i.e., when asymp=TRUE):
#EH = expected value based on sample size
#VH = variance (adjusted for ties if necessary)
#p.value = asymptotic p value
#Sometimes 

[R] Customizing breaks for histograms of unequal ranges

2006-06-30 Thread e . rapsomaniki
Hi, 

I would really appreciate any suggestions on this (rather trivial) problem..

Say I have two vectors:
v1=seq(1:10)
v2=seq(1:15)

For a set of common breaks I need to divide the density of v2 over v1. This
means that I want to avoid having 0 counts for any v1 breakpoint. But
(unsuprisingly) if  I define my common breaks as those returned by calling hist
for v1:

v1.hist=hist(v1, plot=F)
v2.hist=hist(v2, plot=F, breaks=v1.hist$breaks)

I get an error:
Error in hist.default(v1, plot = F, breaks = v2.hist$breaks) : 
some 'x' not counted; maybe 'breaks' do not span range of 'x'

If I had used the combined vector (c(v1,v2)) to set my breaks I would end up
with 0 for some v1 counts. So my question is: is there a way to define breaks
that cover the whole range of v1 and v2 while avoiding having 0 for the
shortest vector?

Many thanks
Eleni Rapsomaniki
Birkbeck College, UK

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


Re: [R] lme convergence

2006-06-30 Thread Dimitris Rizopoulos
I think the following part of lme.formula

if (numIter  controlvals$maxIter) {
stop(Maximum number of iterations reached without convergence.)
}


should be something like


if (numIter  controlvals$maxIter) {
if (controlvals$returnObject) {
warning(Maximum number of iterations reached without 
convergence.)
break
} else {
stop(Maximum number of iterations reached without 
convergence.)
}
}


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


- Original Message - 
From: Spencer Graves [EMAIL PROTECTED]
To: Ravi Varadhan [EMAIL PROTECTED]
Cc: 'Pryseley Assam' [EMAIL PROTECTED]; 'R-Users' 
R-help@stat.math.ethz.ch; Douglas Bates [EMAIL PROTECTED]
Sent: Friday, June 30, 2006 4:08 AM
Subject: Re: [R] lme convergence


   Does anyone know how to obtain the 'returnObject' from an 'lme' 
 run
 that fails to converge?  An argument of this name is described on 
 the
 'lmeControl' help page as, a logical value indicating whether the
 fitted object should be returned when the maximum number of 
 iterations
 is reached without convergence of the algorithm. Default is 
 'FALSE'.

   Unfortunately, I've so far been unable to get it to work, as
 witnessed by the following modification of an example from the 
 '?lme'
 help page:

  library(nlme)
  fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1))
 Error in lme.formula(distance ~ age, data = Orthodont, control =
 lmeControl(msMaxIter = 1)) :
 iteration limit reached without convergence (9)
  fm1
 Error: object fm1 not found
  fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1,
 +  returnObject=TRUE))
 Error in lme.formula(distance ~ age, data = Orthodont, control =
 lmeControl(msMaxIter = 1,  :
 iteration limit reached without convergence (9)
  fm1
 Error: object fm1 not found

   I might be able to fix the problem myself, working through the 
 'lme'
 code line by line, e.g., using 'debug'.  However, I'm not ready to 
 do
 that just now.

   Best Wishes,
   Spencer Graves

 Ravi Varadhan wrote:
 Use try to capture error messages without breaking the loop.
 ?try

 --
 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:r-help-
 [EMAIL PROTECTED] On Behalf Of Pryseley Assam
 Sent: Wednesday, June 28, 2006 12:18 PM
 To: R-Users
 Subject: [R] lme convergence

 Dear R-Users,

   Is it possible to get the covariance matrix from an lme model 
 that did
 not converge ?

   I am doing a simulation which entails fitting linear mixed 
 models, using
 a for loop.
   Within each loop, i generate a new data set and analyze it using 
 a mixed
 model.  The loop stops When the lme function does not converge 
 for a
 simulated dataset. I want to inquire if there is a method to 
 suppress the
 error message from the lme function, or better still, a way of 
 going about
 this issue of the loop ending once the lme function does not 
 converge.

   Thanks in advance,
   Pryseley


 -


 [[alternative HTML version deleted]]

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

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

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


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


Re: [R] GARCH message 75

2006-06-30 Thread Spencer Graves
  For those who can find no documention on the @ operator, ?@ 
says it is used to Extract tbe contents of a slot in a object with a 
formal class structure.  ?slot provides more information, including 
some critical references.  The primary reference is Chambers (1998) 
Programming with Data (Springer).  This describes the S 4 standard for 
object oriented programming in R.  Additional information can be found 
in a chapter of Venables and Ripley (2000) S Programming (Springer).

  When I read Chambers a few years ago, I found it very difficult, 
partly because of the inconsistencies between what Chambers described 
and what was actually implemented in versions of S-Plus and to which I 
had access at the time.  That may have improved since then.

  Hope this helps,
  Spencer Graves

Joe Byers wrote:
 All,
 
 You might check out the @ operator on Classes and objects.  I can find 
 no documentation on this but if you look at code in fseries or fbasic in 
 the method fARMA.summary.  You will find where the object fit in the 
 returned object is access using the @ operator.
 
object - [EMAIL PROTECTED]
ans$coefmat - cbind(format(object$coef,digits=digits), 
 format(object$se.coef,digits=digits),
format(tval,digits=digits), 
 prob=format.pval(prob,digits=digits))
 
 where x is the from x-armaFit().  This I believe would be the same 
 for the GARCH results.
 
 Note the format.pval on the prob.  This is important because if you do 
 not do this you will get 0 for small pvalues.  I learned this by looking 
 at the code of the pretty printing functions so I can save only results 
 that I want from multiple models runs.  You also can then to wald test 
 and LL ratio tests on the models.
 
 Good luck
 Joe
 
 
 
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

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


Re: [R] Running R

2006-06-30 Thread Pramod Anugu
Does any one have solution for R not running? I am getting the error
message. Please Advice.
Thanks

 
   1. gunzip  R-patched.tar.gz
   2. tar -xvf R-patched.tar
   3. changed the directory to the newly created directory R-patched
   4. Typed ./configure
   5. Typed make
   6. Make check
   7. make check-all
   8. Typed make install
   *9. Typed R
# R
Fatal error: unable to open the base package

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


[R] [R-pkgs] twang - Toolkit for Weighting and Analysis of Nonequivalent Groups

2006-06-30 Thread Ridgeway, Greg

The Toolkit for Weighting and Analysis of Nonequivalent Groups (twang
1.0) has been released to CRAN. The package collects functions useful
for computing propensity score weights for treatment effect estimation,
developing nonresponse weights, and diagnosing the quality of those
weights. The package includes a vignette containing some basic theory
and walks through two examples. It is available by typing
vignette(twang)
at the R prompt after loading the library.

Background on the methodology with an application to evaluating drug
treatment programs is available in

McCaffrey, D., G. Ridgeway, A. Morral (2004). Propensity score
estimation with boosted regression for evaluating adolescent substance
abuse treatment, Psychological Methods 9(4):403-425.




This email message is for the sole use of the intended recip...{{dropped}}

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

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


[R] [R-pkgs] gmp: bigintegers with matrix computation

2006-06-30 Thread Antoine
Dear all,

gmp package is now available on cran with version 0.4

Aim of gmp is to provide a lot a function to manipulate big integers, big 
rationals and last but not least big integers in Z/nZ (modulo n)

We add in last version support for matrix computation (standards operators, + * 
- / %*% ...) and inversion matrix (solve):
in Z/nZ if matrix defined in Z/nZ
with rationals (absolute precision) if not.

Best Regards,

Antoine Lucas.

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

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


[R] incomplete final line found by readLines on ...

2006-06-30 Thread Taka Matzmoto
Dear R-users

I need to read some text files produced by some other software.
I used readLines (with n = -1 ) command and then tried to find some numbers 
I liked to extract or some line numbers I like to know.

The problem is that there is no empty line at the end of the text files.

R gave me this warning below;

In addition: Warning message:
incomplete final line found by readLines on 'output.txt'

I know this doesn't affect anything and just warning messages. Is there any 
way to prevent this warning message.

Thanks

Taka

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


[R] Query : Chi Square goodness of fit test

2006-06-30 Thread priti desai
I want to calculate chi square test of goodness of fit to test,
Sample coming from Poisson distribution.

please copy this script in R  run the script
The R script is as follows

## start
#

No_of_Frouds-
c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7)


N - length(No_of_Frouds)

# Estimation of Parameter


lambda- sum(No_of_Frouds)/N
lambda 

pmf  - dpois(i, lambda, log = FALSE)

step_function  - ppois(i, lambda, lower.tail = TRUE, log.p = FALSE)

# Chi-Squared Goodness of Fit Test

# Ho: The data follow a Poisson distribution Vs H1: Not Ho


Frauds - c(1:13)

counts-  c(2,3,3,5,7,2,1,1,2,3,2,1,1,0)  # Observed frequency

Expected
-c(0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.7817821
03,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.2344006
79,0.095299266,0.035764993)

chisq.test(counts, Expected, simulate.p.value =FALSE, correct = FALSE)



# end 


The result of R is as follows

  Pearson's Chi-squared test

data:  counts and poisson_fit 
X-squared = 70, df = 65, p-value = 0.3135

Warning message:
Chi-squared approximation may be incorrect in: chisq.test(counts,
poisson_fit, simulate.p.value = FALSE, correct = FALSE)



But I have done calculations in Excel. I am getting different answer.

Observed  = 2,3,3,5,7,2,1,1,2,3,2,1,1,0
Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78
1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23
4400679,0.095299266,0.035764993


 Estimated Parameter  =4.878788

Chi square stat =  0.000113


My excel answer tally with the book which I have refer for excel.   
Please tell me the correct calculations in R.

##

Awaiting your positive reply.

Regards.
Priti.

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


Re: [R] Running R

2006-06-30 Thread Liaw, Andy
Have you checked and see if there is any error in steps 5 through 8?  If you
can't start R, I doubt steps 6 and 7 ran fine.

Andy 

From: Pramod Anugu 
 
 Does any one have solution for R not running? I am getting 
 the error message. Please Advice.
 Thanks
 
  
    1. gunzip  R-patched.tar.gz
2. tar -xvf R-patched.tar
    3. changed the directory to the newly created directory R-patched
    4. Typed ./configure
    5. Typed make
    6. Make check
    7. make check-all
    8. Typed make install
    *9. Typed R
 # R
 Fatal error: unable to open the base package
 


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


Re: [R] Running R

2006-06-30 Thread Dirk Eddelbuettel
On Fri, Jun 30, 2006 at 07:56:41AM -0500, Pramod Anugu wrote:
 Does any one have solution for R not running? I am getting the error
 message. Please Advice.

Yes, you have emailing the group fairly regularly about your troubles.

And you have been told just about each and every time to *check the output
of the configure step*.  Unless that step concludes with all relevant pieces
found, *there is no point in resuming with steps 5 to 9*.

As you are emailing from an educational institution, consider getting local
help about configuring and compiling software.  Building R is *extensively
documented in a dedicated manual* but the manual may need someone with more
familiarity with the process than you currently have.

Good luck, and consider to stop sending the *same* email to the list every
couple of days.

Dirk


 Thanks
 
 ?
  ? 1. gunzip  R-patched.tar.gz
2. tar -xvf R-patched.tar
 ?? 3. changed the directory to the newly created directory R-patched
 ?? 4. Typed ./configure
 ? ?5. Typed make
 ?? 6. Make check
 ?? 7. make check-all
 ?? 8. Typed make install
 ?? *9. Typed R
 # R
 Fatal error: unable to open the base package
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Hell, there are no rules here - we're trying to accomplish something. 
  -- Thomas A. Edison

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


[R] Random numbers from noncentral t-distribution

2006-06-30 Thread Long Qu
Hi there: 
   
  I'd thought these two versions of noncentral t-distribution are essentially 
the same: 
   qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3))

  But, the scales of the x-axis and the y-axis are quite different according to 
the QQ-plot. 
   
  Did I make any mistakes somewhere? 
   
   
  Thanks, 
  Long
   


-
 Mp3·è¿ñËÑ-иèÈȸè¸ßËÙÏ   
[[alternative HTML version deleted]]

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

Re: [R] incomplete final line found by readLines on ...

2006-06-30 Thread François Pinard
[Taka Matzmoto]

Is there any way to prevent [this] warning message.

Hi, Taka.  The easiest might be using the suppressWarnings wrapper.
See ?suppressWarnings for more information.

-- 
François Pinard   http://pinard.progiciels-bpi.ca

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


Re: [R] Query : Chi Square goodness of fit test

2006-06-30 Thread Jacques VESLOT
  chisq.test(counts, p=Expected/sum(Expected), simulate.p.value =FALSE, 
  correct = FALSE)

 Chi-squared test for given probabilities

data:  counts
X-squared = 40.5207, df = 13, p-value = 0.0001139

Warning message:
l'approximation du Chi-2 est peut-être incorrecte in: chisq.test(counts, p = 
Expected/sum(Expected), 
simulate.p.value = FALSE,

but the use of Chi2 test is incorrect since some of Expected frequencies are 
lower than 5.

---
Jacques VESLOT

CNRS UMR 8090
I.B.L (2ème étage)
1 rue du Professeur Calmette
B.P. 245
59019 Lille Cedex

Tel : 33 (0)3.20.87.10.44
Fax : 33 (0)3.20.87.10.31

http://www-good.ibl.fr
---


priti desai a écrit :
 I want to calculate chi square test of goodness of fit to test,
 Sample coming from Poisson distribution.
 
 please copy this script in R  run the script
 The R script is as follows
 
 ## start
 #
 
 No_of_Frouds-
 c(4,1,6,9,9,10,2,4,8,2,3,0,1,2,3,1,3,4,5,4,4,4,9,5,4,3,11,8,12,3,10,0,7)
 
 
 N - length(No_of_Frouds)
 
 # Estimation of Parameter
 
 
 lambda- sum(No_of_Frouds)/N
 lambda 
 
 pmf  - dpois(i, lambda, log = FALSE)
 
 step_function  - ppois(i, lambda, lower.tail = TRUE, log.p = FALSE)
 
 # Chi-Squared Goodness of Fit Test
 
 # Ho: The data follow a Poisson distribution Vs H1: Not Ho
 
 
 Frauds - c(1:13)
 
 counts-  c(2,3,3,5,7,2,1,1,2,3,2,1,1,0)  # Observed frequency
 
 Expected
 -c(0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.7817821
 03,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.2344006
 79,0.095299266,0.035764993)
 
 chisq.test(counts, Expected, simulate.p.value =FALSE, correct = FALSE)
 
 
 
 # end 
 
 
 The result of R is as follows
 
   Pearson's Chi-squared test
 
 data:  counts and poisson_fit 
 X-squared = 70, df = 65, p-value = 0.3135
 
 Warning message:
 Chi-squared approximation may be incorrect in: chisq.test(counts,
 poisson_fit, simulate.p.value = FALSE, correct = FALSE)
 
 
 
 But I have done calculations in Excel. I am getting different answer.
 
 Observed  = 2,3,3,5,7,2,1,1,2,3,2,1,1,0
 Expected=0.251005528,1.224602726,2.987288468,4.85811559,5.925428863,5.78
 1782103,4.701348074,3.276697142,1.998288788,1.083247457,0.528493456,0.23
 4400679,0.095299266,0.035764993
 
 
  Estimated Parameter  =4.878788
 
 Chi square stat =  0.000113
 
 
 My excel answer tally with the book which I have refer for excel.   
 Please tell me the correct calculations in R.
 
 ##
 
 Awaiting your positive reply.
 
 Regards.
 Priti.
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


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


[R] Way to convert data frame to matrix

2006-06-30 Thread Wade Wall
  I have a text file that I have imported into R.  It contains 3 columns and
316940 rows.  The first column is vegetation plot ID, the second species
names and the third is a cover value (numeric).  I imported using the
read.table function.

My problem is this.  I need to reformat the information as a matrix, with
the first column becoming the row labels and the second the column labels
and the cover values as the matrix cell data.  However, since the
read.tablefunction imported the data as an indexed data frame, I can't
use the columns
as vectors.  Is there a way around this, to convert the data frame as 3
separate vectors?  I have been looking all over for a function, and my
programming skills are not great.

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


Re: [R] Way to convert data frame to matrix

2006-06-30 Thread Mike Sears
?as.matrix

On Friday 30 June 2006 09:11, Wade Wall wrote:
   I have a text file that I have imported into R.  It contains 3 columns
 and 316940 rows.  The first column is vegetation plot ID, the second
 species names and the third is a cover value (numeric).  I imported using
 the read.table function.

 My problem is this.  I need to reformat the information as a matrix, with
 the first column becoming the row labels and the second the column labels
 and the cover values as the matrix cell data.  However, since the
 read.tablefunction imported the data as an indexed data frame, I can't
 use the columns
 as vectors.  Is there a way around this, to convert the data frame as 3
 separate vectors?  I have been looking all over for a function, and my
 programming skills are not great.

 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

-- 
Michael W. Sears, Ph.D.
Associate Editor, Herpetologica
Assistant Professor
Center for Ecology
Department of Zoology
Southern Illinois University
Carbondale, IL 62901

phone: 618-453-4137
web:http://equinox.unr.edu/homepage/msears
http://www.ecology.siu.edu

Natural selection is a mechanism for generating an exceedingly 
high degree of improbability  Sir Ronald A. Fisher (1890-1962)

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


Re: [R] Way to convert data frame to matrix

2006-06-30 Thread Duncan Murdoch
On 6/30/2006 10:11 AM, Wade Wall wrote:
   I have a text file that I have imported into R.  It contains 3 columns and
 316940 rows.  The first column is vegetation plot ID, the second species
 names and the third is a cover value (numeric).  I imported using the
 read.table function.
 
 My problem is this.  I need to reformat the information as a matrix, with
 the first column becoming the row labels and the second the column labels
 and the cover values as the matrix cell data.  However, since the
 read.tablefunction imported the data as an indexed data frame, I can't
 use the columns
 as vectors.  Is there a way around this, to convert the data frame as 3
 separate vectors?  I have been looking all over for a function, and my
 programming skills are not great.

Internally, dataframes are just lists with a class=dataframe 
attribute.  This means you can extract the columns as if they were just 
lists.

So if your columns are named A, B, and C, and the dataframe is dataf, 
you get them as vectors using

dataf$A, dataf$B, and dataf$C

Duncan Murdoch

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


Re: [R] Way to convert data frame to matrix

2006-06-30 Thread Philipp Pagel
On Fri, Jun 30, 2006 at 10:11:15AM -0400, Wade Wall wrote:
   I have a text file that I have imported into R.  It contains 3 columns and
 316940 rows.  The first column is vegetation plot ID, the second species
 names and the third is a cover value (numeric).  I imported using the
 read.table function.
 
 My problem is this.  I need to reformat the information as a matrix, with
 the first column becoming the row labels and the second the column labels
 and the cover values as the matrix cell data.

I'm not 100% sure but I think you are looking for reshape().

cu
Philipp

-- 
Dr. Philipp PagelTel.  +49-8161-71 2131
Dept. of Genome Oriented Bioinformatics  Fax.  +49-8161-71 2186
Technical University of Munich
Science Center Weihenstephan
85350 Freising, Germany

 and

Institute for Bioinformatics / MIPS  Tel.  +49-89-3187 3675
GSF - National Research Center   Fax.  +49-89-3187 3585
  for Environment and Health
Ingolstädter Landstrasse 1
85764 Neuherberg, Germany
http://mips.gsf.de/staff/pagel

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


Re: [R] Way to convert data frame to matrix

2006-06-30 Thread Sarah Goslee
Hi Wade,

On 6/30/06, Wade Wall [EMAIL PROTECTED] wrote:
   I have a text file that I have imported into R.  It contains 3 columns and
 316940 rows.  The first column is vegetation plot ID, the second species
 names and the third is a cover value (numeric).  I imported using the
 read.table function.

 My problem is this.  I need to reformat the information as a matrix, with
 the first column becoming the row labels and the second the column labels
 and the cover values as the matrix cell data.

Try crosstab(mydata$plotID, mydata$species, mydata$cover) from the
ecodist package.

Sarah
-- 
Sarah Goslee

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


[R] Passing arguments to glm()

2006-06-30 Thread Christian Bieli
Hi there

I want to pass arguments (i.e. the response variable and the subset 
argument) in a self-made function to glm.
Here is one way I can do this:

f.myglm - function(y,subfact,subval) {
  
glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval)
}

  str(d.mydata)
`data.frame':15806 obs. of  3 variables:
 $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ...
 $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ...
 $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ...

  f.myglm('y','x2','yes')

But is there a way I can pass the arguments and use the data argument of 
glm()?
In a naive way of thinking I'd like to something like this:
f.myglm - function(y,sub) {
  glm(y~x1,family=binomial,data=d.mydata,subset=sub)
}
  f.myglm(y=y,sub=x2=='yes')

I know that's not possible, because the objects y and x2 are not defined 
in the user workspace.
So, something like passing the arguments as an expression and evaluate 
it in the glm function should work, but I didn't manage to do it.

I'd appreciate your advice.
Christian

  R.version
 _ 
platform i386-pc-mingw32
arch i386  
os   mingw32   
system   i386, mingw32 
status 
major2 
minor2.1   
year 2005  
month12
day  20
svn rev  36812 
language R

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


[R] apply a function to several lists' components

2006-06-30 Thread Taka Matzmoto
Dear R-user
I have 100 lists.
Each list has several components.
For example,

data1
$a
[1] 1 2

$b
[1] 3 4

$c
[1] 5

There are data1, data2,, data100. All lists have the same number and the 
same name of components.


Is there any function I can use for applying to only a specific component 
across 100 lists?
(e.g.,  taking mean of $c acorss 100 lists) or do I need to write my own 
function for that?

Thank you.

Taka,

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


Re: [R] Passing arguments to glm()

2006-06-30 Thread Jacques VESLOT
  f.myglm - function(y=y, subset=x2 == 'yes', data=d.d.mydata) 
  eval(parse(text=glm(, 
deparse(substitute(y)), ~ x1, family=binomial, data=, 
deparse(substitute(data)), , subset =, 
subset, )))
  f.myglm()

---
Jacques VESLOT

CNRS UMR 8090
I.B.L (2ème étage)
1 rue du Professeur Calmette
B.P. 245
59019 Lille Cedex

Tel : 33 (0)3.20.87.10.44
Fax : 33 (0)3.20.87.10.31

http://www-good.ibl.fr
---


Christian Bieli a écrit :
 Hi there
 
 I want to pass arguments (i.e. the response variable and the subset 
 argument) in a self-made function to glm.
 Here is one way I can do this:
 
 f.myglm - function(y,subfact,subval) {
   
 glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval)
 }
 
   str(d.mydata)
 `data.frame':15806 obs. of  3 variables:
  $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ...
  $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ...
  $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ...
 
   f.myglm('y','x2','yes')
 
 But is there a way I can pass the arguments and use the data argument of 
 glm()?
 In a naive way of thinking I'd like to something like this:
 f.myglm - function(y,sub) {
   glm(y~x1,family=binomial,data=d.mydata,subset=sub)
 }
   f.myglm(y=y,sub=x2=='yes')
 
 I know that's not possible, because the objects y and x2 are not defined 
 in the user workspace.
 So, something like passing the arguments as an expression and evaluate 
 it in the glm function should work, but I didn't manage to do it.
 
 I'd appreciate your advice.
 Christian
 
   R.version
  _ 
 platform i386-pc-mingw32
 arch i386  
 os   mingw32   
 system   i386, mingw32 
 status 
 major2 
 minor2.1   
 year 2005  
 month12
 day  20
 svn rev  36812 
 language R
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


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


Re: [R] Passing arguments to glm()

2006-06-30 Thread Jacques VESLOT
i forgot a paste():
f.myglm - function(y=y, subset=x2 == 'yes', data=d.d.mydata) 
eval(parse(text=paste(glm(,
deparse(substitute(y)), ~ x1, family=binomial, data=, 
deparse(substitute(data)), , subset =,
subset, ), sep=)))
---
Jacques VESLOT

CNRS UMR 8090
I.B.L (2ème étage)
1 rue du Professeur Calmette
B.P. 245
59019 Lille Cedex

Tel : 33 (0)3.20.87.10.44
Fax : 33 (0)3.20.87.10.31

http://www-good.ibl.fr
---

Jacques VESLOT a écrit :
   f.myglm - function(y=y, subset=x2 == 'yes', data=d.d.mydata) 
 eval(parse(text=glm(, 
 deparse(substitute(y)), ~ x1, family=binomial, data=, 
 deparse(substitute(data)), , subset =, 
 subset, )))
   f.myglm()
 
 ---
 Jacques VESLOT
 
 CNRS UMR 8090
 I.B.L (2ème étage)
 1 rue du Professeur Calmette
 B.P. 245
 59019 Lille Cedex
 
 Tel : 33 (0)3.20.87.10.44
 Fax : 33 (0)3.20.87.10.31
 
 http://www-good.ibl.fr
 ---
 
 
 Christian Bieli a écrit :
 
Hi there

I want to pass arguments (i.e. the response variable and the subset 
argument) in a self-made function to glm.
Here is one way I can do this:

f.myglm - function(y,subfact,subval) {
  
glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval)
}

  str(d.mydata)
`data.frame':15806 obs. of  3 variables:
 $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ...
 $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ...
 $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ...

  f.myglm('y','x2','yes')

But is there a way I can pass the arguments and use the data argument of 
glm()?
In a naive way of thinking I'd like to something like this:
f.myglm - function(y,sub) {
  glm(y~x1,family=binomial,data=d.mydata,subset=sub)
}
  f.myglm(y=y,sub=x2=='yes')

I know that's not possible, because the objects y and x2 are not defined 
in the user workspace.
So, something like passing the arguments as an expression and evaluate 
it in the glm function should work, but I didn't manage to do it.

I'd appreciate your advice.
Christian

  R.version
 _ 
platform i386-pc-mingw32
arch i386  
os   mingw32   
system   i386, mingw32 
status 
major2 
minor2.1   
year 2005  
month12
day  20
svn rev  36812 
language R

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

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


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


[R] median of gamma distribution

2006-06-30 Thread Philip He
Doese anyone know a R function to find the median of a gamma distribution?

[[alternative HTML version deleted]]

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


Re: [R] median of gamma distribution

2006-06-30 Thread Peter Dalgaard
Philip He [EMAIL PROTECTED] writes:

 Doese anyone know a R function to find the median of a gamma distribution?

qgamma(0.5, )

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

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


Re: [R] median of gamma distribution

2006-06-30 Thread Thomas Lumley
On Fri, 30 Jun 2006, Philip He wrote:

 Doese anyone know a R function to find the median of a gamma distribution?


It's not clear what you mean.  If you know the parameters of a gamma 
distribution then qgamma() will give you any quantile.

If you have data and want to estimate the median then it's hard to beat 
median(), but you could use mle() to estimate the parameters and then 
qgamma().

-thomas

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


Re: [R] median of gamma distribution

2006-06-30 Thread markleeds
From: Philip He [EMAIL PROTECTED]
Date: Fri Jun 30 10:30:28 CDT 2006
To: R help list r-help@stat.math.ethz.ch
Subject: [R] median of gamma distribution


someone might know a  more elegant way but one way 
is to use the distribution ( i think it's dnorm for
the normal but i get confused because ther
is also pnorm, rnorm and one other but you can look that up to
figure out which one is the one you need )
functions in R and just put in .50 as the probabiliy
input and itr will kick back the result which is the median.

  mark

   


 

Doese anyone know a R function to find the median of a gamma distribution?

   [[alternative HTML version deleted]]

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

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


Re: [R] median of gamma distribution

2006-06-30 Thread Ted Harding
On 30-Jun-06 Philip He wrote:
 Doese anyone know a R function to find the median of a gamma
 distribution?

qgamma will do it. Test:

 -log(0.5)
[1] 0.6931472
 qgamma(0.5,1)
[1] 0.6931472

Ted.



E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 30-Jun-06   Time: 16:53:16
-- XFMail --

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


Re: [R] getting the smoother matrix from smooth.spline

2006-06-30 Thread Simon Wood
smooth.matrix = function(x, df){
  n = length(x);
  A = matrix(0, n, n);
  for(i in 1:n){
y = rep(0, n); y[i]=1;
yi = predict(smooth.spline(x, y, df=df),x)$y;
A[,i]= yi;
}
  (A+t(A))/2;
}


- Simon Wood, Mathematical Sciences, University of Bath, Bath BA2 7AY
- +44 (0)1225 386603 www.maths.bath.ac.uk/~sw283/


On Sat, 24 Jun 2006, Gregory Gentlemen wrote:

 Can anyone tell me the trick for obtaining the smoother matrix from 
 smooth.spline when there are non-unique values for x. I have the following 
 code but, of course, it only works when all values of x are unique.

  ## get the smoother matrix (x having unique values
 smooth.matrix = function(x, df){
 n = length(x);
 A = matrix(0, n, n);
 for(i in 1:n){
   y = rep(0, n); y[i]=1;
   yi = smooth.spline(x, y, df=df)$y;
   A[,i]= yi;
 }
 (A+t(A))/2;
 }


  Thanks for any assistance,
  Gregory


 -

 -
 Get a sneak peak at messages with a handy reading pane.
   [[alternative HTML version deleted]]

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


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


Re: [R] Random numbers from noncentral t-distribution

2006-06-30 Thread Thomas Lumley
On Fri, 30 Jun 2006, Long Qu wrote:

 Hi there:

  I'd thought these two versions of noncentral t-distribution are essentially 
 the same:
   qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3))

  But, the scales of the x-axis and the y-axis are quite different according 
 to the QQ-plot.

  Did I make any mistakes somewhere?


No, I think we did.

We have
 rt
function (n, df, ncp = 0)
{
 if (ncp == 0)
 .Internal(rt(n, df))
 else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df))
}

and the rchisq() in the denominator should be inside the sqrt().


-thomas

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


Re: [R] getting the smoother matrix from smooth.spline

2006-06-30 Thread Gabor Grothendieck
Perhaps this could be developed into a spline smooth method
for model.matrix and included in R.

On 6/30/06, Simon Wood [EMAIL PROTECTED] wrote:
 smooth.matrix = function(x, df){
  n = length(x);
  A = matrix(0, n, n);
  for(i in 1:n){
y = rep(0, n); y[i]=1;
yi = predict(smooth.spline(x, y, df=df),x)$y;
A[,i]= yi;
 }
  (A+t(A))/2;
 }


 - Simon Wood, Mathematical Sciences, University of Bath, Bath BA2 7AY
 - +44 (0)1225 386603 www.maths.bath.ac.uk/~sw283/


 On Sat, 24 Jun 2006, Gregory Gentlemen wrote:

  Can anyone tell me the trick for obtaining the smoother matrix from 
  smooth.spline when there are non-unique values for x. I have the following 
  code but, of course, it only works when all values of x are unique.
 
   ## get the smoother matrix (x having unique values
  smooth.matrix = function(x, df){
  n = length(x);
  A = matrix(0, n, n);
  for(i in 1:n){
y = rep(0, n); y[i]=1;
yi = smooth.spline(x, y, df=df)$y;
A[,i]= yi;
  }
  (A+t(A))/2;
  }
 
 
   Thanks for any assistance,
   Gregory
 
 
  -
 
  -
  Get a sneak peak at messages with a handy reading pane.
[[alternative HTML version deleted]]
 
  __
  R-help@stat.math.ethz.ch mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide! 
  http://www.R-project.org/posting-guide.html
 

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


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


Re: [R] Passing arguments to glm()

2006-06-30 Thread Gabor Grothendieck
This is the same as glm except that
- formula may be of class character in which case its
  regarded as the name of the response variable and
  the formula defaults to resp ~ Species for that response
- the data frame defaults to iris
- modify as appropriate for your case

myglm - function (formula, family = gaussian, data,
weights, subset, na.action, start = NULL, etastart, mustart, offset,
control = glm.control(...), model = TRUE, method = glm.fit,
x = FALSE, y = TRUE, contrasts = NULL, ...) {
cl - match.call()
cl[[1]] - as.name(glm)
if (is.character(formula)) {
   fo - . ~ Species  ### default formula
   fo[[2]] - as.name(formula)
   cl$formula - fo
}
if (missing(data)) cl$data - as.name(iris) # default data frame
eval(cl, parent.frame())
}
# test
myglm(Sepal.Length, subset = Petal.Length  mean(Petal.Length))
myglm(Sepal.Length ~ Petal.Length)
myglm(Sepal.Length, subset = 1:100)


On 6/30/06, Christian Bieli [EMAIL PROTECTED] wrote:
 Hi there

 I want to pass arguments (i.e. the response variable and the subset
 argument) in a self-made function to glm.
 Here is one way I can do this:

 f.myglm - function(y,subfact,subval) {

 glm(d.mydata[,y]~d.mydata[,'x1'],family=binomial,subset=d.mydata[,subfact]==subval)
 }

   str(d.mydata)
 `data.frame':15806 obs. of  3 variables:
  $ y : Factor w/ 2 levels no,yes: 1 1 1 1 1 1 1 NA 1 1 ...
  $ x1: Factor w/ 2 levels no,yes: 2 2 1 2 2 2 2 2 2 2 ...
  $ x2: Factor w/ 2 levels no,yes: 1 1 1 1 1 2 2 1 2 2 ...

   f.myglm('y','x2','yes')

 But is there a way I can pass the arguments and use the data argument of
 glm()?
 In a naive way of thinking I'd like to something like this:
 f.myglm - function(y,sub) {
  glm(y~x1,family=binomial,data=d.mydata,subset=sub)
 }
   f.myglm(y=y,sub=x2=='yes')

 I know that's not possible, because the objects y and x2 are not defined
 in the user workspace.
 So, something like passing the arguments as an expression and evaluate
 it in the glm function should work, but I didn't manage to do it.

 I'd appreciate your advice.
 Christian

   R.version
 _
 platform i386-pc-mingw32
 arch i386
 os   mingw32
 system   i386, mingw32
 status
 major2
 minor2.1
 year 2005
 month12
 day  20
 svn rev  36812
 language R

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


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


Re: [R] lme convergence

2006-06-30 Thread Michael Jerosch-Herold
It looks like in the call to lme

 fm1 - lme(distance ~ age, data = Orthodont,
+control=lmeControl(msMaxIter=1))

you did not specify any random effects. Why not try:

 fm1 - lme(distance ~ age, random= ~1| groupID, data = Orthodont,
+control=lmeControl(msMaxIter=1))

where groupID is some factor that can be used to stratify the data.

Also, the Othodont data set is used in Pinheiro  Bates book, and you may 
want to consult
that book to see the models they use in connection with that data set. For the 
Orthodont data set
the groupID would most likely be the subject ID (Subject variable).

So a possible model would be:

 fm1 - lme(distance ~ age, random= ~1|Subject, data=Orthodont)
 summary(fm1)
Linear mixed-effects model fit by REML
 Data: Orthodont 
   AIC  BIClogLik
  455.0025 465.6563 -223.5013

Random effects:
 Formula: ~1 | Subject
(Intercept) Residual
StdDev:2.114724 1.431592

Fixed effects: distance ~ age 
Value Std.Error DF  t-value p-value
(Intercept) 16.76 0.8023952 80 20.5   0
age  0.660185 0.0616059 80 10.71626   0
 Correlation: 
(Intr)
age -0.845

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max 
-3.66453932 -0.53507984 -0.01289591  0.48742859  3.72178465 

Number of Observations: 108
Number of Groups: 27 

So this runs fine.

As, I said this data set and its analysis is discussed extensively in Pinheiro 
and Bates book

Michael Jerosch-Herold


 Spencer Graves [EMAIL PROTECTED] 06/29/06 7:08 PM 
  Does anyone know how to obtain the 'returnObject' from an 'lme' run 
that fails to converge?  An argument of this name is described on the 
'lmeControl' help page as, a logical value indicating whether the 
fitted object should be returned when the maximum number of iterations 
is reached without convergence of the algorithm. Default is 'FALSE'.

  Unfortunately, I've so far been unable to get it to work, as 
witnessed by the following modification of an example from the '?lme' 
help page:

  library(nlme)
  fm1 - lme(distance ~ age, data = Orthodont,
+control=lmeControl(msMaxIter=1))
Error in lme.formula(distance ~ age, data = Orthodont, control = 
lmeControl(msMaxIter = 1)) :
iteration limit reached without convergence (9)
  fm1
Error: object fm1 not found
  fm1 - lme(distance ~ age, data = Orthodont,
+control=lmeControl(msMaxIter=1,
+  returnObject=TRUE))
Error in lme.formula(distance ~ age, data = Orthodont, control = 
lmeControl(msMaxIter = 1,  :
iteration limit reached without convergence (9)
  fm1
Error: object fm1 not found   

  I might be able to fix the problem myself, working through the 'lme' 
code line by line, e.g., using 'debug'.  However, I'm not ready to do 
that just now.

  Best Wishes,
  Spencer Graves

Ravi Varadhan wrote:
 Use try to capture error messages without breaking the loop.
 ?try
 
 --
 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:r-help-
 [EMAIL PROTECTED] On Behalf Of Pryseley Assam
 Sent: Wednesday, June 28, 2006 12:18 PM
 To: R-Users
 Subject: [R] lme convergence

 Dear R-Users,

   Is it possible to get the covariance matrix from an lme model that did
 not converge ?

   I am doing a simulation which entails fitting linear mixed models, using
 a for loop.
   Within each loop, i generate a new data set and analyze it using a mixed
 model.  The loop stops When the lme function does not converge for a
 simulated dataset. I want to inquire if there is a method to suppress the
 error message from the lme function, or better still, a way of going about
 this issue of the loop ending once the lme function does not converge.

   Thanks in advance,
   Pryseley


 -


  [[alternative HTML version deleted]]

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

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


Re: [R] Random numbers from noncentral t-distribution

2006-06-30 Thread Long Qu
Thank you very much for your kind reply. It solved the problem of rt( ). :D
   
  But it seems that the qt( ) also have problems: 
  I modified the rt( ) function as you suggested,
  rt - function (n, df, ncp = 0)
{
if (ncp == 0)
.Internal(rt(n, df))
else rnorm(n, ncp)/sqrt(rchisq(n, df)/df)
}

  Then I increase the number of random variables to 1, and made a QQ-plot:
   qqplot(rt(1,df=20,ncp=3),qt(runif(1),df=20,ncp=3))

  I've got some spurious points at lower-left corner. It seems that qt( ) 
results were truncated. 
   
  I also tried this with another df and ncp: 
   pt(-.75,df=2,ncp=1)
[1] 0.05726429
 sum(qt(1:1/10001,df=2,ncp=1)  -.75)
[1] 0
where I'd expected the last number should be  550 or so, not 0. 
  
 
   
  Thanks again, the modified rt( ) is now OK for my work. 
  Long 
   
  
Thomas Lumley [EMAIL PROTECTED] wrote£º
  On Fri, 30 Jun 2006, Long Qu wrote:

 Hi there:

 I'd thought these two versions of noncentral t-distribution are essentially 
 the same:
  qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3))

 But, the scales of the x-axis and the y-axis are quite different according to 
 the QQ-plot.

 Did I make any mistakes somewhere?


No, I think we did.

We have
 rt
function (n, df, ncp = 0)
{
if (ncp == 0)
.Internal(rt(n, df))
else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df))
}

and the rchisq() in the denominator should be inside the sqrt().


-thomas



-
ÇÀ×¢ÑÅ»¢Ãâ·ÑÓÊÏä-3.5GÈÝÁ¿£¬20M¸½¼þ£¡ 
[[alternative HTML version deleted]]

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

Re: [R] lme convergence

2006-06-30 Thread Doran, Harold
In the old version of lme, one could construct a grouped data object and
this would alleviate the need to specify the random portion of the
model. So, Spencer's call is equivalent to 

fm1 - lme(distance ~ age, random= ~age| Subject, data = Orthodont)

This condition does not hold under lmer, however.

 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of 
 Michael Jerosch-Herold
 Sent: Friday, June 30, 2006 12:37 PM
 To: [EMAIL PROTECTED]
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] lme convergence
 
 It looks like in the call to lme
 
  fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1))
 
 you did not specify any random effects. Why not try:
 
  fm1 - lme(distance ~ age, random= ~1| groupID, data = Orthodont,
 +control=lmeControl(msMaxIter=1))
 
 where groupID is some factor that can be used to stratify the data.
 
 Also, the Othodont data set is used in Pinheiro  Bates 
 book, and you may want to consult that book to see the models 
 they use in connection with that data set. For the Orthodont 
 data set the groupID would most likely be the subject ID 
 (Subject variable).
 
 So a possible model would be:
 
  fm1 - lme(distance ~ age, random= ~1|Subject, data=Orthodont)
  summary(fm1)
 Linear mixed-effects model fit by REML
  Data: Orthodont 
AIC  BIClogLik
   455.0025 465.6563 -223.5013
 
 Random effects:
  Formula: ~1 | Subject
 (Intercept) Residual
 StdDev:2.114724 1.431592
 
 Fixed effects: distance ~ age 
 Value Std.Error DF  t-value p-value
 (Intercept) 16.76 0.8023952 80 20.5   0
 age  0.660185 0.0616059 80 10.71626   0
  Correlation: 
 (Intr)
 age -0.845
 
 Standardized Within-Group Residuals:
 Min  Q1 Med  Q3 Max 
 -3.66453932 -0.53507984 -0.01289591  0.48742859  3.72178465 
 
 Number of Observations: 108
 Number of Groups: 27 
 
 So this runs fine.
 
 As, I said this data set and its analysis is discussed 
 extensively in Pinheiro and Bates book
 
 Michael Jerosch-Herold
 
 
  Spencer Graves [EMAIL PROTECTED] 06/29/06 7:08 PM 
 Does anyone know how to obtain the 'returnObject' 
 from an 'lme' run that fails to converge?  An argument of 
 this name is described on the 'lmeControl' help page as, a 
 logical value indicating whether the fitted object should be 
 returned when the maximum number of iterations is reached 
 without convergence of the algorithm. Default is 'FALSE'.
 
 Unfortunately, I've so far been unable to get it to 
 work, as witnessed by the following modification of an 
 example from the '?lme' 
 help page:
 
   library(nlme)
   fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1))
 Error in lme.formula(distance ~ age, data = Orthodont, 
 control = lmeControl(msMaxIter = 1)) :
   iteration limit reached without convergence (9)   fm1
 Error: object fm1 not found
   fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1,
 +  returnObject=TRUE))
 Error in lme.formula(distance ~ age, data = Orthodont, 
 control = lmeControl(msMaxIter = 1,  :
   iteration limit reached without convergence (9)   fm1
 Error: object fm1 not found 
 
 I might be able to fix the problem myself, working 
 through the 'lme' 
 code line by line, e.g., using 'debug'.  However, I'm not 
 ready to do that just now.
 
 Best Wishes,
 Spencer Graves
 
 Ravi Varadhan wrote:
  Use try to capture error messages without breaking the loop.
  ?try
  
  
 --
  
  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:r-help- 
  [EMAIL PROTECTED] On Behalf Of Pryseley Assam
  Sent: Wednesday, June 28, 2006 12:18 PM
  To: R-Users
  Subject: [R] lme convergence
 
  Dear R-Users,
 
Is it possible to get the covariance matrix from an lme 
 model that 
  did not converge ?
 
I am doing a simulation which entails fitting linear 
 mixed models, 
  using a for loop.
Within each loop, i generate a new data set and analyze 
 it using a 
  mixed model.  The loop stops When the lme function does not 
  converge for a simulated dataset. I want to inquire if there is a 
  method to suppress the error message from the lme 
 function, or better 
  still, a way of going about this issue of the loop ending 
 once the lme function does not converge.
 
Thanks in advance,
Pryseley
 
 
  -
 
 
 

Re: [R] Random numbers from noncentral t-distribution

2006-06-30 Thread Thomas Lumley
On Sat, 1 Jul 2006, Long Qu wrote:

 Thank you very much for your kind reply. It solved the problem of rt( ). :D

  But it seems that the qt( ) also have problems:

Yes, there does seem to be a problem near zero. A clearer version is
   curve(qt(x,df=20,ncp=3),from=0,to=0.004)
   curve(qt(10^x,df=20,ncp=3),from=-10,to=-2,n=1000)

The fix is less obvious here. I'll file it as a bug.

-thomas

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

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


Re: [R] Running R

2006-06-30 Thread Liaw, Andy
Just for the record.  (One of the) problem was that configure picked up
ATLAS, but had problem with it at link time for whatever reason.

(This is on some version of Redhat, x86_64.  I should think there are people
who have similar setup and got both ATLAS and readline to work.)

Andy

 From: Pramod Anugu
 
 IT WORKS!!!
 Thanks you very much for your help. So the below line fixed 
 my installation
 
 #./configure --with-blas=no --with-readline=no
 


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


[R] tkbutton command - how to know which button was clicked?

2006-06-30 Thread JeeBee

In the below code fragment, print(arg) always prints the
last element of rekeningen$rekening.
Is this because of lazy evaluation? I.e. arg is evaluated at
the time the button is pressed?
And, if so, how can I avoid this?
I tried function() {force(arg); print(arg)} but that didn't work either.

Thanks,
Jeebee.

  for(rek in seq(1,nrow(rekeningen))) {
arg - rekeningen$rekening[rek]

tkgrid(tkbutton(frame.1,
  text=paste(Saldo historie, arg),
  command=function() print(arg)),
  sticky=news)
  }

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


Re: [R] lme convergence

2006-06-30 Thread Spencer Graves
  Harold is correct.  The help page for 'Orthodont' includes the 
following example:

  formula(Orthodont)
distance ~ age | Subject

  If 'random' is not specified, 'lme' sets random = formula(data).  If 
that's NULL, the 'lme' help page says it Defaults to a formula 
consisting of the right hand side of 'fixed'.  This will generally 
return an error, indicated by the following example:

  tstDF - data.frame(gp=rep(1:2, 2), y=1:4)
  lme(y~1, tstDF)
Error in getGroups.data.frame(dataMix, groups) :
Invalid formula for groups
  lme(y~1, tstDF, random=~1)
Error in getGroups.data.frame(dataMix, groups) :
Invalid formula for groups
 
  Hope this helps.
  Spencer

Doran, Harold wrote:
 In the old version of lme, one could construct a grouped data object and
 this would alleviate the need to specify the random portion of the
 model. So, Spencer's call is equivalent to 
 
 fm1 - lme(distance ~ age, random= ~age| Subject, data = Orthodont)
 
 This condition does not hold under lmer, however.
 
 -Original Message-
 From: [EMAIL PROTECTED] 
 [mailto:[EMAIL PROTECTED] On Behalf Of 
 Michael Jerosch-Herold
 Sent: Friday, June 30, 2006 12:37 PM
 To: [EMAIL PROTECTED]
 Cc: r-help@stat.math.ethz.ch
 Subject: Re: [R] lme convergence

 It looks like in the call to lme

  fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1))

 you did not specify any random effects. Why not try:

  fm1 - lme(distance ~ age, random= ~1| groupID, data = Orthodont,
 +control=lmeControl(msMaxIter=1))

 where groupID is some factor that can be used to stratify the data.

 Also, the Othodont data set is used in Pinheiro  Bates 
 book, and you may want to consult that book to see the models 
 they use in connection with that data set. For the Orthodont 
 data set the groupID would most likely be the subject ID 
 (Subject variable).

 So a possible model would be:

 fm1 - lme(distance ~ age, random= ~1|Subject, data=Orthodont)
 summary(fm1)
 Linear mixed-effects model fit by REML
  Data: Orthodont 
AIC  BIClogLik
   455.0025 465.6563 -223.5013

 Random effects:
  Formula: ~1 | Subject
 (Intercept) Residual
 StdDev:2.114724 1.431592

 Fixed effects: distance ~ age 
 Value Std.Error DF  t-value p-value
 (Intercept) 16.76 0.8023952 80 20.5   0
 age  0.660185 0.0616059 80 10.71626   0
  Correlation: 
 (Intr)
 age -0.845

 Standardized Within-Group Residuals:
 Min  Q1 Med  Q3 Max 
 -3.66453932 -0.53507984 -0.01289591  0.48742859  3.72178465 

 Number of Observations: 108
 Number of Groups: 27 

 So this runs fine.

 As, I said this data set and its analysis is discussed 
 extensively in Pinheiro and Bates book

 Michael Jerosch-Herold


 Spencer Graves [EMAIL PROTECTED] 06/29/06 7:08 PM 
Does anyone know how to obtain the 'returnObject' 
 from an 'lme' run that fails to converge?  An argument of 
 this name is described on the 'lmeControl' help page as, a 
 logical value indicating whether the fitted object should be 
 returned when the maximum number of iterations is reached 
 without convergence of the algorithm. Default is 'FALSE'.

Unfortunately, I've so far been unable to get it to 
 work, as witnessed by the following modification of an 
 example from the '?lme' 
 help page:

   library(nlme)
   fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1))
 Error in lme.formula(distance ~ age, data = Orthodont, 
 control = lmeControl(msMaxIter = 1)) :
  iteration limit reached without convergence (9)   fm1
 Error: object fm1 not found
   fm1 - lme(distance ~ age, data = Orthodont,
 +control=lmeControl(msMaxIter=1,
 +  returnObject=TRUE))
 Error in lme.formula(distance ~ age, data = Orthodont, 
 control = lmeControl(msMaxIter = 1,  :
  iteration limit reached without convergence (9)   fm1
 Error: object fm1 not found

I might be able to fix the problem myself, working 
 through the 'lme' 
 code line by line, e.g., using 'debug'.  However, I'm not 
 ready to do that just now.

Best Wishes,
Spencer Graves

 Ravi Varadhan wrote:
 Use try to capture error messages without breaking the loop.
 ?try


 --
 
 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:r-help- 
 [EMAIL PROTECTED] On Behalf Of Pryseley Assam
 Sent: Wednesday, June 28, 2006 12:18 PM
 To: R-Users
 Subject: [R] lme 

Re: [R] apply a function to several lists' components

2006-06-30 Thread JeeBee

Maybe this helps

( data1 = list(a=c(1,2), b=c(3,4), c=c(5,6,7)) )
( data2 = list(a=c(10,11), b=c(30,40), c=c(70,80)) )

cc - NULL
for(data in ls(pattern=^data[0-9]+$)) {
  cc - c(cc, with(get(data), c))
}

mean(cc)

JeeBee.


On Fri, 30 Jun 2006 09:50:51 -0500, Taka Matzmoto wrote:

 Dear R-user
 I have 100 lists.
 Each list has several components.
 For example,
 
data1
 $a
 [1] 1 2
 
 $b
 [1] 3 4
 
 $c
 [1] 5
 
 There are data1, data2,, data100. All lists have the same number and the 
 same name of components.
 
 
 Is there any function I can use for applying to only a specific component 
 across 100 lists?
 (e.g.,  taking mean of $c acorss 100 lists) or do I need to write my own 
 function for that?
 
 Thank you.
 
 Taka,
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

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


[R] send output to printer

2006-06-30 Thread Matthias Braeunig
It has to be a simple thing, but I could not figure it out:

How do I send the text output from object x to the printer?
As a shell user I would expect a pipe to the printer... |kprinter or
|lpr -Pmyprinter somehow. And yes, I'm on Linux.

Thanks!

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


[R] SAS Proc Mixed and lme

2006-06-30 Thread Kellie J. Archer, Ph.D.
I am trying to use lme to fit a mixed effects model to get the same
results as when using the following SAS code:

proc mixed;
class refseqid probeid probeno end;
model expression=end logpgc / ddfm=satterth;
random probeno probeid / subject=refseqid type=cs;
lsmeans end / diff cl; run;

There are 3 genes (refseqid) which is the large grouping factor, with
2 probeids nested within each refseqid, and 16 probenos nested within
each of the probeids.

I have specified in the SAS Proc Mixed procedure that the
variance-covariance structure is to be compound symmetric. Therefore,
the variance-covariance matrix is a block diagonal matrix of the form,

V_1  0   0
0   V_2  0
00   V3

where each V_i represents a RefSeqID. Moreover, for each V_i the
structure within the block is

v_{11}   v{12}
v_{21}   v{22}

where v_{11} and v_{22} are different probeids nested within the
refseqid, and so are correlated. The structure of
both v_{11} and v_{22} are compound symmetric, and v_{12} and v{21}
contain a constant for all elements of the matrix.

I have tried to reproduce this using lme, but it is unclear from the
documentation (and Pinheiro  Bates text) how the pdBlocked and
compound symmetric structure can be combined.

fit.lme-lme(expression~End+logpgc,random=list(RefSeqID=pdBlocked(list
(~1,~ProbeID-1),pdClass=pdSymm)),data=dataset,correlation=corCompSym
m(form=~1|RefSeqID/ProbeID/ProbeNo))


The point estimates are essentially the same comparing R and SAS for
the fixed effects, but the 95% confidence intervals are much shorter
using lme(). In order to find the difference in the algorithms used by
SAS and R I tried to extract the variance-covariance matrix to look at
its structure. I used the getVarCov() command, but it tells me that
this function is not available for nested structures. Is there another
way to extract the variance-covariance structure for nested models?
Does anyone know how I could get the var-cov structure above using
lme?


Kellie J. Archer, Ph.D.
Assistant Professor, Department of Biostatistics
Fellow, Center for the Study of Biological Complexity
Virginia Commonwealth University
1101 East Marshall St., B1-066
Richmond, VA 23298-0032
phone: (804) 827-2039
fax: (804) 828-8900
e-mail: [EMAIL PROTECTED]
website: www.people.vcu.edu/~kjarcher

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


[R] lme and SAS Proc mixed

2006-06-30 Thread Kellie J. Archer, Ph.D.
I am trying to use lme to fit a mixed effects model to get the same
results as when using the following SAS code:

proc mixed;
class refseqid probeid probeno end;
model expression=end logpgc / ddfm=satterth;
random probeno probeid / subject=refseqid type=cs;
lsmeans end / diff cl; run;

There are 3 genes (refseqid) which is the large grouping factor, with
2 probeids nested within each refseqid, and 16 probenos nested within
each of the probeids.

I have specified in the SAS Proc Mixed procedure that the
variance-covariance structure is to be compound symmetric. Therefore,
the variance-covariance matrix is a block diagonal matrix of the form,

V_1  0   0
0   V_2  0
00   V3

where each V_i represents a RefSeqID. Moreover, for each V_i the
structure within the block is 

v_{11}   v{12}
v_{21}   v{22}

where v_{11} and v_{22} are different probeids nested within the
refseqid, and so are correlated. The structure of
both v_{11} and v_{22} are compound symmetric, and v_{12} and v{21}
contain a constant for all elements of the matrix.

I have tried to reproduce this using lme, but it is unclear from the
documentation (and Pinheiro  Bates text) how the pdBlocked and
compound symmetric structure can be combined.

fit.lme-lme(expression~End+logpgc,random=list(RefSeqID=pdBlocked(list
(~1,~ProbeID-1),pdClass=pdSymm)),data=dataset,correlation=corCompSym
m(form=~1|RefSeqID/ProbeID/ProbeNo))


The point estimates are essentially the same comparing R and SAS for
the fixed effects, but the 95% confidence intervals are much shorter
using lme(). In order to find the difference in the algorithms used by
SAS and R I tried to extract the variance-covariance matrix to look at
its structure. I used the getVarCov() command, but it tells me that
this function is not available for nested structures. Is there another
way to extract the variance-covariance structure for nested models?
Does anyone know how I could get the var-cov structure above using
lme?


Kellie J. Archer, Ph.D.
Assistant Professor, Department of Biostatistics
Fellow, Center for the Study of Biological Complexity
Virginia Commonwealth University
1101 East Marshall St., B1-066
Richmond, VA 23298-0032
phone: (804) 827-2039
fax: (804) 828-8900
e-mail: [EMAIL PROTECTED]
website: www.people.vcu.edu/~kjarcher

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


[R] Creating Vectors

2006-06-30 Thread Raphael Fraser
 type count
   0 20
   1 15
0 10
   1 35
   0 28

I would like to create two vectors from the data above. For example,
type1=c(15, 35) and type0 = c(20, 10, 28). Can any one help

Raphael

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


Re: [R] tkbutton command - how to know which button was clicked?

2006-06-30 Thread Peter Dalgaard
JeeBee [EMAIL PROTECTED] writes:

 In the below code fragment, print(arg) always prints the
 last element of rekeningen$rekening.
 Is this because of lazy evaluation? I.e. arg is evaluated at
 the time the button is pressed?

No and yes. Lazy evaluation has nothing to do with it, but the
function passed to command= is not evaluated until the button is
pressed. If you want a different arg variable for each
button-function, you have to make sure that they have different
environments or have the arg actually part of the function.

There are various possible variations. One is

mk.button - function(arg)
 tkgrid(tkbutton(frame.1,
   text=paste(Saldo historie, arg),
   command=function() print(arg)),
   sticky=news)

for  {

mk.button(arg)
}

another is to wrap a local({arg - arg;.}) around the tkgrid call
in your code. A 3rd one is

.command=eval(substitute(function()print(arg)))...

 And, if so, how can I avoid this?
 I tried function() {force(arg); print(arg)} but that didn't work either.
 
 Thanks,
 Jeebee.
 
   for(rek in seq(1,nrow(rekeningen))) {
 arg - rekeningen$rekening[rek]
 
 tkgrid(tkbutton(frame.1,
   text=paste(Saldo historie, arg),
   command=function() print(arg)),
   sticky=news)
   }

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


[R] useR! 2006: presentation slides

2006-06-30 Thread Achim Zeileis
Dear useRs,

the useR! 2006 conference took place in Vienna two weeks ago: it was an
exciting and interesting meeting with about 400 useRs from all over the
world and more than 150 presentations.

Especially for those of you who could not make it to the conference, we
have made 4up PDF versions of the presentations slides available.

The slides for the keynote lectures are available at
  http://www.R-project.org/useR-2006/Keynotes/
and the user-contributed presentations at
  http://www.R-project.org/useR-2006/Presentations/

Finally, some materials for the panel discussion on Getting recognition
for excellence in computational statistics are provided at
  http://www.R-project.org/useR-2006/PanelDisc/
including a summary, the short presentation slides of the panelists and
the full discussion as a video.

A big thank you to everyone who contributed to the conference and...
best wishes from Vienna!

The organizing team
Achim, Torsten, David, Bettina, Fritz and Kurt.

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

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


Re: [R] Creating Vectors

2006-06-30 Thread Alexander Nervedi
Hi all.

I have a factor variable distributed over time. I am looking for an elegant 
way to code duration of a state. Suppose,

rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, prob 
= unit.p),
+  label = c(Drought, Normal, High))

rainfall.shocks
[1] Normal  HighHighDrought Normal  Normal  HighNormal  Drought
[10] Normal  Drought Normal  Normal  Normal  Normal


So capture the duration of say drought, I'd need a variable that is able to 
keep track of rainfall.shocks as well as its past values. I was wondering if 
there is any obvious way to do this. the Drought variable in this case would 
have values

0 0 0 1 0 0 0 0 1 0 1 0 0 0 0

many thanks for the suggestions you are likely to make.

Alexander Nervedi

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


[R] Fwd: time series patterns

2006-06-30 Thread Raphael Fraser
-- Forwarded message --
From: Alexander Nervedi [EMAIL PROTECTED]
Date: Jun 30, 2006 4:56 PM
Subject: time series patterns
To: [EMAIL PROTECTED]


Hi all.

I have a factor variable distributed over time. I am looking for an elegant
way to code duration of a state. Suppose,

rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, prob
= unit.p),
+  label = c(Drought, Normal, High))

rainfall.shocks
[1] Normal  HighHighDrought Normal  Normal  HighNormal  Drought
[10] Normal  Drought Normal  Normal  Normal  Normal


So capture the duration of say drought, I'd need a variable that is able to
keep track of rainfall.shocks as well as its past values. I was wondering if
there is any obvious way to do this. the Drought variable in this case would
have values

0 0 0 1 0 0 0 0 1 0 1 0 0 0 0

many thanks for the suggestions you are likely to make.

Alexander Nervedi

_
Express yourself instantly with MSN Messenger! Download today - it's FREE!
http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/

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


Re: [R] {Spam?} {Virus?} Message could not be delivered

2006-06-30 Thread J-Francois Paillotin
I am out of the office and will be back on Monday 31 July.
 
Your mail will not be forwarded automatically.
Your mail regarding {Spam?} {Virus?} Message could not be delivered will be 
read when I return

In urgent cases please contact Mr Kholdoun Torki

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


[R] weird error message

2006-06-30 Thread Alexander Nervedi
Hi!

In the example below why is x[[1]] == 0.2237724 false?

Alexander Nervedi


x - runif(10)
x[[1]]
[1] 0.2237724

x
[1] 0.2237724 0.2678944 0.9375811 0.5963889 0.6180519 0.6449580 0.7308510
[8] 0.7347386 0.4837286 0.1416100

x[[1]] == 0.2237724
FALSE






From: Alexander Nervedi [EMAIL PROTECTED]
To: r-help@stat.math.ethz.ch
Subject: Re: [R] Creating Vectors
Date: Fri, 30 Jun 2006 21:57:35 +

Hi all.

I have a factor variable distributed over time. I am looking for an elegant
way to code duration of a state. Suppose,

 rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE, 
prob
 = unit.p),
+  label = c(Drought, Normal, High))
 
 rainfall.shocks
[1] Normal  HighHighDrought Normal  Normal  HighNormal  Drought
[10] Normal  Drought Normal  Normal  Normal  Normal


So capture the duration of say drought, I'd need a variable that is able to
keep track of rainfall.shocks as well as its past values. I was wondering 
if
there is any obvious way to do this. the Drought variable in this case 
would
have values

0 0 0 1 0 0 0 0 1 0 1 0 0 0 0

many thanks for the suggestions you are likely to make.

Alexander Nervedi

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

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


Re: [R] Creating Vectors

2006-06-30 Thread jim holtman
Does this do what you want?  It creates a 'list' with the vectors:


 x - 'type count
+   0 20
+   1 15
+0 10
+   1 35
+   0 28
+ '
 x - read.table(textConnection(x), header=TRUE)
 x
  type count
1020
2115
3010
4135
5028
 type - split(x$count, x$type)
 type
$`0`
[1] 20 10 28

$`1`
[1] 15 35

 type[['0']]
[1] 20 10 28
 type[['1']]
[1] 15 35




On 6/30/06, Raphael Fraser [EMAIL PROTECTED] wrote:

 type count
   0 20
   1 15
0 10
   1 35
   0 28

 I would like to create two vectors from the data above. For example,
 type1=c(15, 35) and type0 = c(20, 10, 28). Can any one help

 Raphael

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




-- 
Jim Holtman
Cincinnati, OH
+1 513 646 9390 (Cell)
+1 513 247 0281 (Home)

What is the problem you are trying to solve?

[[alternative HTML version deleted]]

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


Re: [R] weird error message

2006-06-30 Thread Gabor Grothendieck
This is FAQ 7.31:

http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f

Also please do not piggy back on other threads since it makes the
archives less useful.


On 6/30/06, Alexander Nervedi [EMAIL PROTECTED] wrote:
 Hi!

 In the example below why is x[[1]] == 0.2237724 false?

 Alexander Nervedi


 x - runif(10)
 x[[1]]
 [1] 0.2237724

 x
 [1] 0.2237724 0.2678944 0.9375811 0.5963889 0.6180519 0.6449580 0.7308510
 [8] 0.7347386 0.4837286 0.1416100

 x[[1]] == 0.2237724
 FALSE






 From: Alexander Nervedi [EMAIL PROTECTED]
 To: r-help@stat.math.ethz.ch
 Subject: Re: [R] Creating Vectors
 Date: Fri, 30 Jun 2006 21:57:35 +
 
 Hi all.
 
 I have a factor variable distributed over time. I am looking for an elegant
 way to code duration of a state. Suppose,
 
  rainfall.shocks - factor(sample(c(1,2,3), size = 15, replace = TRUE,
 prob
  = unit.p),
 +  label = c(Drought, Normal, High))
  
  rainfall.shocks
 [1] Normal  HighHighDrought Normal  Normal  HighNormal  Drought
 [10] Normal  Drought Normal  Normal  Normal  Normal
 
 
 So capture the duration of say drought, I'd need a variable that is able to
 keep track of rainfall.shocks as well as its past values. I was wondering
 if
 there is any obvious way to do this. the Drought variable in this case
 would
 have values
 
 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0
 
 many thanks for the suggestions you are likely to make.
 
 Alexander Nervedi
 
 __
 R-help@stat.math.ethz.ch mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide!
 http://www.R-project.org/posting-guide.html

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


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


[R] postscript file too large : maybe an R question

2006-06-30 Thread markleeds
i created a postscipt file in R and then i downloaded a free version
of ghostview to view it. unfortunately, i get the message

fata error : dynamic memory exhausted
when i try to view it.

when i do a dir on windows xp, the file size is 149,034,475
and i know there about 17,000 graphs. is there
a way of possibly viewing this size postscript file in R itself ?

   Thanks

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


[R] polynomial expansion in R

2006-06-30 Thread yyan liu
Hi:
   I have two vectors of data, x and y and I want to get the polynomial 
expansion of (x+y)^p with any integer power p in R. Suppose p=2, then I want a 
matrix of five vectors, namely, x y x^2 y^2 x*y. The coefficient of the 
polynomial is not needed. I can write it manully if p is small. But I want it 
in the case of p=10 or even bigger, is there any function in R can do that 
automatically for me with certain choice of p?
   Thx a lot!
 
 liu
 

-


[[alternative HTML version deleted]]

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


[R] generate bi-variate normal data

2006-06-30 Thread Shin, David
Dear all,

I would like to generate bi-variate normal data given that the first column
of the data is known. for example:
I first generate a set of data using the command, 
x - rmvnorm(10, c(0, 0), matrix(c(1, 0, 0, 1), 2))

then I would like to sum up the two columns of x:
x.sum - apply(x, 1, sum)

now with x.sum I would like to generate another column of data, say y, that
makes cbind(x.sum, y) follow a bi-variate normal distribution with mean =
c(0, 0) and sigma = matrix(c(1, 0, 0, 1),2)

I will appreciate for all insights.

David s.

 
This email may contain confidential material.\ If you were n...{{dropped}}

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


Re: [R] generate bi-variate normal data

2006-06-30 Thread markleeds
From: Shin, David [EMAIL PROTECTED]
Date: Sat Jul 01 00:15:21 CDT 2006
To: 'r-help@stat.math.ethz.ch' r-help@stat.math.ethz.ch
Subject: [R] generate bi-variate normal data

it's an interesting question. someone else
on this list can answer more explicitly but
i think you have to use the result for the multivariate
normal distribution ( bivariate case ) where , if the
joint is normal , then the conditional is normal also
with parameters a function of the 2 means and the elements of
the covariance matrix. the result in any decent mathematical statistics such as 
 casella berger. so, given the one column, generate the other column 
conditionally using the formula  and then the joint dist will be bivariate 
normal.




Dear all,

I would like to generate bi-variate normal data given that the first column
of the data is known. for example:
I first generate a set of data using the command, 
x - rmvnorm(10, c(0, 0), matrix(c(1, 0, 0, 1), 2))

then I would like to sum up the two columns of x:
x.sum - apply(x, 1, sum)

now with x.sum I would like to generate another column of data, say y, that
makes cbind(x.sum, y) follow a bi-variate normal distribution with mean =
c(0, 0) and sigma = matrix(c(1, 0, 0, 1),2)

I will appreciate for all insights.

David s.

 
This email may contain confidential material.\ If you were n...{{dropped}}

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

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