Re: [R] main title x title and y title with ggplot2

2008-03-05 Thread ONKELINX, Thierry
Guillaume,

Have a look at the ggplot book on p. 29
(http://had.co.nz/ggplot2/book.pdf). 

HTH,

Thierry



ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
[EMAIL PROTECTED] 
www.inbo.be 

Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney

-Oorspronkelijk bericht-
Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Namens guillaume chaumet
Verzonden: woensdag 5 maart 2008 8:51
Aan: r-help@r-project.org
Onderwerp: [R] main title x title and y title with ggplot2

Hi R people,
I'm a R newbie and I'm trying to put main title, x title and y title in
my
graph with no success.
Any idea?
I'm sorry for this newbie question..

Guillaume

[[alternative HTML version deleted]]

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

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


[R] Marginal Effects in a Logit Regression

2008-03-05 Thread Paul Sweeting
Hi

 

I am carrying out some logit regressions, so have a (0,1) dependent variable
and various dependent variables, each in the (-inf, inf) range.  As well as
the usual output, I need to report the marginal effects of each the
explanatory variables.  How can this be done?  I've seen a couple of similar
questions on the mailing list but no answers.

 

Thanks!

 

Paul


[[alternative HTML version deleted]]

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


Re: [R] main title x title and y title with ggplot2

2008-03-05 Thread guillaume chaumet
Thierry
First thank you for the celerity of your response.
Second I use ggplot2 like this :
ggplot(data, aes(x,y,fill)) + geom_point() + etc.
Where did you your xlab and ylab when using ggplot2 like that?

Guillaume

2008/3/5, ONKELINX, Thierry [EMAIL PROTECTED]:

 Guillaume,

 Have a look at the ggplot book on p. 29
 (http://had.co.nz/ggplot2/book.pdf).

 HTH,

 Thierry

 
 
 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek / Research Institute for Nature
 and Forest
 Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
 methodology and quality assurance
 Gaverstraat 4
 9500 Geraardsbergen
 Belgium
 tel. + 32 54/436 185
 [EMAIL PROTECTED]
 www.inbo.be

 Do not put your faith in what statistics say until you have carefully
 considered what they do not say.  ~William W. Watt
 A statistical analysis, properly conducted, is a delicate dissection of
 uncertainties, a surgery of suppositions. ~M.J.Moroney

 -Oorspronkelijk bericht-
 Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
 Namens guillaume chaumet
 Verzonden: woensdag 5 maart 2008 8:51
 Aan: r-help@r-project.org
 Onderwerp: [R] main title x title and y title with ggplot2


 Hi R people,
 I'm a R newbie and I'm trying to put main title, x title and y title in
 my
 graph with no success.
 Any idea?
 I'm sorry for this newbie question..

 Guillaume


 [[alternative HTML version deleted]]

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


[[alternative HTML version deleted]]

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


Re: [R] main title x title and y title with ggplot2

2008-03-05 Thread ONKELINX, Thierry
Guillaume,

You'll have to add the appropriate scales. 

ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your
xlabel) + scale_y_continuous(your ylabel)

I suppose you can add a main title in a similar way, but I haven't found
that yet. But I shure that Hadley will answer this.

Thierry



ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
[EMAIL PROTECTED] 
www.inbo.be 

Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney

-Oorspronkelijk bericht-
Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Namens guillaume chaumet
Verzonden: woensdag 5 maart 2008 9:45
Aan: r-help@r-project.org
Onderwerp: Re: [R] main title x title and y title with ggplot2

Thierry
First thank you for the celerity of your response.
Second I use ggplot2 like this :
ggplot(data, aes(x,y,fill)) + geom_point() + etc.
Where did you your xlab and ylab when using ggplot2 like that?

Guillaume

2008/3/5, ONKELINX, Thierry [EMAIL PROTECTED]:

 Guillaume,

 Have a look at the ggplot book on p. 29
 (http://had.co.nz/ggplot2/book.pdf).

 HTH,

 Thierry



 
 ir. Thierry Onkelinx
 Instituut voor natuur- en bosonderzoek / Research Institute for Nature
 and Forest
 Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
 methodology and quality assurance
 Gaverstraat 4
 9500 Geraardsbergen
 Belgium
 tel. + 32 54/436 185
 [EMAIL PROTECTED]
 www.inbo.be

 Do not put your faith in what statistics say until you have carefully
 considered what they do not say.  ~William W. Watt
 A statistical analysis, properly conducted, is a delicate dissection
of
 uncertainties, a surgery of suppositions. ~M.J.Moroney

 -Oorspronkelijk bericht-
 Van: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED]
 Namens guillaume chaumet
 Verzonden: woensdag 5 maart 2008 8:51
 Aan: r-help@r-project.org
 Onderwerp: [R] main title x title and y title with ggplot2


 Hi R people,
 I'm a R newbie and I'm trying to put main title, x title and y title
in
 my
 graph with no success.
 Any idea?
 I'm sorry for this newbie question..

 Guillaume


 [[alternative HTML version deleted]]

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


[[alternative HTML version deleted]]

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

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


Re: [R] vertex labels in igraph from adjacency matrix

2008-03-05 Thread Gabor Csardi
On Wed, Mar 05, 2008 at 02:27:21AM -0500, Charilaos Skiadas wrote:
[...]
 
 Btw, you will likely want to take the betweenness call out, and call  
 it once and store the result, instead of calling it twice (well,  
 assuming the graph is largish). Or even better, use which.max:
 
 which.max(betweenness(graph = my.graph, v=V(my.graph), directed =  
 FALSE))

This is almost good, but there is a catch, in igraph vertices are 
numbered from zero. So if you want an igraph vertex id, then you 
need to subtract one from this, i.e.:

maxb - which.max(betweennness(my.graph, directed=FALSE))-1

You can double check it:

betweenness(my.graph, maxb, directed=FALSE)

Gabor

PS. there is also an igraph mailing list, see the igraph homepage
at igraph.sf.net

 Haris Skiadas
 Department of Mathematics and Computer Science
 Hanover College
 
[...]

-- 
Csardi Gabor [EMAIL PROTECTED]UNIL DGM

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


[R] ipf function in R

2008-03-05 Thread Chandra Shah
Hi
I have a 3 x 2 contingency table:
10 20
30 40
50 60
I want to update the frequencies to new marginal totals:
100 130
40 80 110
I want to use  the ipf (iterative proportional fitting) function which 
is apparently in the cat package.
Can somebody please advice me how to input this data and invoke ipf in R 
to obtain an updated contingency table?
Thanks.
By the way I am quite new to R.

-- 
 

Dr Chandra Shah
Senior Research Fellow
Monash University-ACER Centre for the Economics of Education and Training
Faculty of Education, Building 6,
Monash University
Victoria
Australia 3800
Tel. +61 3 9905 2787
Fax +61 3 9905 9184

www.education.monash.edu.au/centres/ceet

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


[R] coxme - fitting random treatment effect nested within centre

2008-03-05 Thread Rasanga Ruwanthi
Dear all,

I am using coxme function in Kinship library to fit random treatment effect 
nested within centre. I got 3 treatments (0,1,2) and 3 centres. I used 
following commands, but got an error. 

 ugroup=paste(rep(1:3,each=3),rep(0:2,3),sep='/')
 mat1=bdsmatrix(rep(c(1,1,1,1,1,1,1,1,1),3),blocksize=rep(3,3),dimnames=list(ugroup,ugroup))
 mat2=bdsmatrix(rep(c(0,0,0,0,0,0,0,0,1),3),blocksize=rep(3,3),dimnames=list(ugroup,ugroup))
 group=paste(dat1$centre,dat1$treat,sep='/')
 coxme(Surv(time,status) ~ as.factor(treat), data=dat1,random= 
 ~1|group,varlist=list(mat1,mat2),rescale=F,pdcheck=FALSE)

Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, : 
Random effects variance is not spd
 
Could anyone help me correcting this error? 
Many thanks in advance.
Ruwanthi


  

Looking for last minute shopping deals?

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


Re: [R] vertex labels in igraph from adjacency matrix

2008-03-05 Thread Gabor Csardi
Mark, graph.adjacency always preserves the order of the vertices,
so the vertex at row/column 1 will be vertex #0 in the igraph graph,
etc. I'll document this in a minute. 

This means that you can always do

g - graph.adjacency(A)
V(g)$name - colnames(A)

But i completely agree that this should be done by default,
i'll do that in the near future.

Btw, weighted shortest path calculation (= where the edges have 
weights or capacities) is only implemented in the development 
version of igraph. Just in case you're looking for that.

Best,
Gabor

On Wed, Mar 05, 2008 at 01:39:41AM -0500, Mark W Kimpel wrote:
 I am getting some unexpected results from some functions of igraph and 
 it is possible that I am misinterpreting the vertex numbers. Eg., the 
 max betweenness measure seems to be from a vertex that is not connected 
 to a single other vertex. Below if my code snippet:
 
 require(igraph)
 my.graph - graph.adjacency(adjmatrix = my.adj.matrix, mode=c(undirected))
 
 most.between.vert - colnames(my.adj.matrix)[which(betweenness(graph = 
 my.graph, v=V(my.graph), directed = FALSE) == max(betweenness(graph = 
 my.graph, v=V(my.graph), directed = FALSE))) + 1]
 
 sum(my.adj.matrix[most.between.vertext,])
 0
 
 
 Is there a way to automatically set the vertex name attributes from the 
 row and colnames of the inputted adjacency matrix to graph.adjacency? 
 This would seem like an intuitive feature, but I can't find it 
 documented. In the absence of this feature, how can I make sure that I 
 am setting the vertex name attributes correctly?
 -- 
 
 Mark W. Kimpel MD  ** Neuroinformatics ** Dept. of Psychiatry
 Indiana University School of Medicine
 
 15032 Hunter Court, Westfield, IN  46074
 
 (317) 490-5129 Work,  Mobile  VoiceMail
 (317) 204-4202 Home (no voice mail please)
 
 mwkimpelatgmaildotcom
 
 **

-- 
Csardi Gabor [EMAIL PROTECTED]UNIL DGM

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


Re: [R] main title x title and y title with ggplot2

2008-03-05 Thread Ingo Michaelis


ONKELINX, Thierry Thierry.ONKELINX at inbo.be writes:

 
 Guillaume,
 
 You'll have to add the appropriate scales. 
 
 ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your
 xlabel) + scale_y_continuous(your ylabel)
 
 I suppose you can add a main title in a similar way, but I haven't found
 that yet. But I shure that Hadley will answer this.
 
 Thierry

Guillaume,

I had the same problem and found a solution in some forums. Try this:

p-ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your
xlabel) + scale_y_continuous(your ylabel)

p$title-your title

print(p)

Ingo.

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


[R] questions about p.adjust

2008-03-05 Thread Ng Stanley
Hi,

I have some questions about p.adjust.

The false discovery rate is a less stringent condition than the family wise
error rate, so these methods are more powerful than the others., these
methods refer to FDR methods or FWER methods. Simply what are the
differences/pros/cons of both classes of methods.

The 'BH' and 'BY' methods control the false discovery rate, are they FDR
method ?

fdr stands for false discovery rate, how different is it from BH or BY ?

Lastly, which method should I use concerning microarray analysis for
selecting differentially expressed genes.

TIA

[[alternative HTML version deleted]]

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


[R] queries about bitmap()

2008-03-05 Thread Ng Stanley
Hi,

There are different types of tiff methods in bitmap(), which one should be
used for publication-quality pictures ? 'tiffcrle',
'tiffg3', 'tiffg32d', 'tiffg4', 'tifflzw', 'tiffpack',
'tiff12nc',
 'tiff24nc',

Thanks

[[alternative HTML version deleted]]

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


[R] R_alloc with structures with flexible array members

2008-03-05 Thread Ramon Diaz-Uriarte
Dear All,

In a package, I want to use some C code where I am using a structure
(as the basic element of a linked list) with flexible array members.
Basically, this is a structure where the last component is an
incomplete array type  (e.g., Harbison  Steel, C, a reference
manual, 5th ed., p. 159) such as:

struct Sequence {
  struct Sequence *next;
  int len;
  unsigned int state_count[];
};


To create one such sequence, I allocate storage (following Harbison
and Steel) in a C program as follows:

struct Sequence *A;
int n = 4;
A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int));


If I understand correctly, however, it would be better to use R_alloc
instead of malloc (memory automagically freed on exit and error;
error-checking). But I do not know how to make the call to R_alloc
here, since R_alloc allocates n units of size bytes each.


I've tried, without success, the following two:

int to_add_for_R_alloc =
(int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int));

  A = (struct sequence *) R_alloc(to_add_for_R_alloc + n,
sizeof(unsigned int));

or even a brute force attempt as:

 A = (struct sequence *) R_alloc( 100,  sizeof(struct sequence));


but both result in segmentation faults.


Should I just keep using malloc (and free at end)?

Thanks,

R.
-- 
Ramon Diaz-Uriarte
Statistical Computing Team
Structural Biology and Biocomputing Programme
Spanish National Cancer Centre (CNIO)
http://ligarto.org/rdiaz

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


Re: [R] ipf function in R

2008-03-05 Thread Ted Harding
On 05-Mar-08 07:14:28, Chandra Shah wrote:
 Hi
 I have a 3 x 2 contingency table:
 10 20
 30 40
 50 60
 I want to update the frequencies to new marginal totals:
 100 130
 40 80 110
 I want to use  the ipf (iterative proportional fitting) function
 which is apparently in the cat package.
 Can somebody please advice me how to input this data and invoke
 ipf in R to obtain an updated contingency table?
 Thanks.
 By the way I am quite new to R.

It is not clear from your description what you really want to do.
Hence, not clear whether you should be using ipf() to do it!

The purpose of ipf() in the 'cat' package is very simple:
You suppply a complete contingency table, and specify the
order of interactions you want to include in the fit (using
the margins parameter in ipf()); and then the value of ipf()
is a table of the fitted (expected) values corresponding to
a model of the type you have specified, fitted by maximum
likelihood.

Your task seems to be different. If I understand your statement.
it looks as though you want to find a,b,c,d such that

  10+a  20+c | 40
  30+b  40+d | 80
  5060   |110
-+---
 100   130   |230

or, since there are row and column constraints on a,b,c,d:

  10+a  30-a | 40
  40-a  40+a | 80
  5060   |110
-+---
 100   130   |230

so then the task would be to find a.

But now it depends on what you would mean by find a.

You could treat 'a' as an unknown parameter, and fit it
by maximum likelihood under various assumptions about the
dependence between rows and columns (e.g. that the log
odds-ratios have particular values).

Or perhaps (and here is where ipf() may come in) for a
given choice of 'a' you could submit the resulting table
to ipf() and compute chisq from the resulting fitted
values and the given observed values; then seek the
value of 'a' which minimises chisq.

However, it's not really possible to advise in more detail
without being sure of what you are trying to do. Above,
I am only guessing!

Best wishes,
Ted.


E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 05-Mar-08   Time: 10:26:57
-- XFMail --

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


Re: [R] queries about bitmap()

2008-03-05 Thread Uwe Ligges


Ng Stanley wrote:
 Hi,
 
 There are different types of tiff methods in bitmap(), which one should be
 used for publication-quality pictures ? 

I'd avoid bitmap at all and try vector formats such as PostScript or PDF 
or similar formats  depending on your publication.

Uwe Ligges


 'tiffcrle',
 'tiffg3', 'tiffg32d', 'tiffg4', 'tifflzw', 'tiffpack',
 'tiff12nc',
  'tiff24nc',
 
 Thanks
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] Marginal Effects in a Logit Regression

2008-03-05 Thread Stefan Grosse
On Wednesday 05 March 2008 09:25:11 am Paul Sweeting wrote:
PS I am carrying out some logit regressions, so have a (0,1) dependent
PS , I need to report the marginal effects of

Have a look at lrm of the design package. It reports Dxy which is maybe what 
you want...

Stefan

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


Re: [R] problem

2008-03-05 Thread Erika Frigo
Goodmorning Jim,
My file has not only more than a million values, but more than a million 
rows and moreless 30 columns (it is a productive dataset for cows), infact 
with read.table i'm not able to import it.
It is an xls file.
How do you import your million rows and 4-5 columns files?
thank you
regards

Dr.ssa Erika Frigo
Università degli Studi di Milano
Facoltà di Medicina Veterinaria
Dipartimento di Scienze e Tecnologie Veterinarie per la Sicurezza Alimentare 
(VSA)

Via Grasselli, 7
20137 Milano
Tel. 02/50318515
Fax 02/50318501
- Original Message - 
From: jim holtman [EMAIL PROTECTED]
To: Erika Frigo [EMAIL PROTECTED]
Cc: r-help@r-project.org
Sent: Tuesday, March 04, 2008 6:13 PM
Subject: Re: [R] problem


 Is it just a file with a million values or is it some type of a
 structure with a million rows of indeterinent columns?  If it is just
 a million numbers, you can easily read with is 'scan' or 'read.table'
 with no problem.  I work with data structures that have several
 million rows and 4-5 columns without any problems.  What is the format
 of the input?

 On 3/4/08, Erika Frigo [EMAIL PROTECTED] wrote:
 Good evening to everybody,
 I have  problems to import in R a really big dataset (more than 100 
 values). Which is the best package to install?
 Is there someone who works with this kind of dataset and can help me, 
 please?

 Thank you very much,
 Regards





 Dr.ssa Erika Frigo
 Department of Veterinary Sciences and Technology for Food Safety
 University of Milan

 Via Grasselli, 7
 20137 Milano
 Tel.  +39 0250318515
 Fax  +39 0250318501
[[alternative HTML version deleted]]

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



 -- 
 Jim Holtman
 Cincinnati, OH
 +1 513 646 9390

 What is the problem you are trying to solve?

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


[R] Principle component analysis

2008-03-05 Thread phthao05

Thanks to Mr.Liviu Androvic and Mr.Richard Rowe helped me in PCA. 
Because I have just learn R language in a few day so I have many problem. 
1) I don't know why PCA rotation function not run although I try many times.
Would you please hepl me and explain how to read the PCA map (both of
rotated and unrotated) in a concrete example. 
2) Where I can find document relate: Plan S(A), S(A*B), S(A)*B?
   Thanks alot.
-- 
View this message in context: 
http://www.nabble.com/Principle-component-analysis-tp15846902p15846902.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Spidergram

2008-03-05 Thread Antony Unwin
A parallel coordinate plot would do fine.  Load the package iplots  
and then use the command ipcp(x1, x2,...)

Antony Unwin

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


Re: [R] R-Terminal

2008-03-05 Thread Martin Kaffanke
Thats really great, but now the sink() uses also this width.  Is it
possible to make sink using another width (i.e. 80 characters)?

Thanks,
Martin


Am Dienstag, den 04.03.2008, 13:23 +0100 schrieb Martin Elff:
 On Tuesday 04 March 2008 (12:34:47), Peter Dalgaard wrote:
  Martin Kaffanke wrote:
   Hi there!
  
   I use an gnome-terminal for using R.  When I resize the termial to the
   maximum size, R uses only the left side of the window.  Can I tell R to
   use the whole window somehow?
 
  This seems to do it:
 
  options(width=Sys.getenv(COLUMNS))
 
  (Surprisingly, at least to me, the COLUMNS environment variable appears
  to be automagically updated by gnome-terminal even after R is started. I
  was expecting to have to play with system(tset) and suchlike.)
 
 That seems to work also with KDE's Konsole. 
 
 It may be worthwile considering
 to implement some automatic width-adjustment like:
 
 .adjustwidth - function(...){
   options(width=Sys.getenv(COLUMNS))
   TRUE
 }
 
 addTaskCallback(.adjustwidth)
 
 After that, the output width is adapted to the actual terminal width
 each time one hits return.
 
 
 Just my 20 cents.
 Best,
 Martin
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
-- 
Ihr Partner für Webdesign, Webapplikationen und Webspace.
http://www.roomandspace.com/
Martin Kaffanke +43 650 4514224


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


[R] Reversed but positive axis in trellis plots?

2008-03-05 Thread Fredrik Karlsson
Hi,

In my discpipline, it is common to plot one acoustic property on a
positive scale but from top to bottom on the ordinate and the same for
another measurement on the abscissa.
So, the origin of the plot is on the top right of the plot, with
increasing values to the left /down. This is to highlight the
correlation between the acoustic measurement and the position of the
forming structure, for instance when teaching it to students.

The grouping ability of the trellis plot is quite handy whan plotting
many instances of the same thing, so I was wondering if it is possible
to make a trellis xyplot behave this way?
Converting all values to negative and changing labels to the negative
of the negative seems one solution to the reverseness of the axes, but
how do I change the position? Is it possible?

/Fredrik

-- 
Give up learning, and put an end to your troubles.

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


[R] non-linear correlation

2008-03-05 Thread Irene Mantzouni
 

Hi all!

 

This is a rather statistical question; 

Which effect sizes (parametric or not) could I use in order to estimate
the amount of non-linear correlation between 2 variables? 

Is it possible to correct for auto-correlation within the correlated
times series? 

Any suggestions for the appropriate packages/ functions are more than
welcome!!

 

Thank you!


[[alternative HTML version deleted]]

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


[R] t.test p-Value

2008-03-05 Thread Eleni Christodoulou
Hello list,

I am trying to apply the paired t.test between diseased and not diseased
patients to identify genes that are more expressed in the one situation
under the other. In order to retrieve the genes that are more expressed in
the positive disease state I do:
p.values-c()
for(i in 1:length(Significant[,1])){
p.values[i]-try(t.test(positive[i,],negative[i,],alternative
=greater)$p.value)
}

which(p.values0.01)


where Significant is my matrix of  genes  and their expression in tumors and
positive, negative are subsets of thes matrix.
Whn p0.01, I reject the null hypothesis and I accept the alternative one,
that I have greater gene expression in positive than in negative.
I assume I must be doing sth wrong because the heatmap that I get with the
genes that pass the filter of p-value is wrong.

Could anyone help me with this?

thanks a lot,
Eleni

[[alternative HTML version deleted]]

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


Re: [R] problem

2008-03-05 Thread Philipp Pagel
On Wed, Mar 05, 2008 at 12:32:19PM +0100, Erika Frigo wrote:
 My file has not only more than a million values, but more than a million 
 rows and moreless 30 columns (it is a productive dataset for cows), infact 
 with read.table i'm not able to import it.
 It is an xls file.

read.table() expects clear text -- e.g. csv or tab separated in the case
of read.delim(). If your file is in xls format the simplest option would
be to export the data to CSV format from Excel.

If for some reason that is not an option please have a look at the R
Data Import/Export manual.

Of course neither will solve the problem of not enough memory if your
file is simply too large. In that case you will may want to put your
data into a database and have R connect to it and retrieve the data in
smaller chunks as required.

cu
Philipp

-- 
Dr. Philipp Pagel  Tel.  +49-8161-71 2131
Lehrstuhl für Genomorientierte Bioinformatik   Fax.  +49-8161-71 2186
Technische Universität München
Wissenschaftszentrum Weihenstephan
85350 Freising, Germany
 
 and
 
Institut für Bioinformatik und Systembiologie / MIPS
Helmholtz Zentrum München -
Deutsches Forschungszentrum für Gesundheit und Umwelt
Ingolstädter Landstrasse 1
85764 Neuherberg, Germany
http://mips.gsf.de/staff/pagel

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


Re: [R] R-Terminal

2008-03-05 Thread Martin Elff
On Wednesday 05 March 2008 (12:56:17), Martin Kaffanke wrote:
 Thats really great, but now the sink() uses also this width.  Is it
 possible to make sink using another width (i.e. 80 characters)?

# auto width adjustment

.adjustWidth - function(...){
   options(width=Sys.getenv(COLUMNS))
   TRUE
}
 
.adjustWidthCallBack - addTaskCallback(.adjustWidth)

# a wrapper for 'sink'

Sink - function(...,width=80){
 if(nargs()){
 removeTaskCallback(.adjustWidthCallBack)
 rm(.adjustWidthCallBack,envir=globalenv())
 options(width=width)
 sink(...)
}
 else{
sink()
assign(.adjustWidthCallBack,
addTaskCallback(.adjustWidth),
envir=globalenv()
)
 }
}

# test
Sink(test.out,width=20) # nonsensically narrow, but shows that
  # width=something has an effect...
iris[1:10,]
Sink()
as.matrix(readLines(test.out))
iris[1:10,]


Best,
Martin


-
Dr. Martin Elff
Faculty of Social Sciences
LSPWIVS (van Deth)
University of Mannheim
A5, 6
68131 Mannheim
Germany

Phone: +49-621-181-2093
Fax: +49-621-181-2099
E-Mail: [EMAIL PROTECTED]
Web: http://webrum.uni-mannheim.de/sowi/elff/
 http://www.sowi.uni-mannheim.de/lspwivs/
-

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


Re: [R] t.test p-Value

2008-03-05 Thread Eleni Christodoulou
I am sorry, the test is unpaired...But my question remains

Thanks,
Eleni

On Wed, Mar 5, 2008 at 2:33 PM, Eleni Christodoulou [EMAIL PROTECTED]
wrote:

 Hello list,

 I am trying to apply the paired t.test between diseased and not diseased
 patients to identify genes that are more expressed in the one situation
 under the other. In order to retrieve the genes that are more expressed in
 the positive disease state I do:
 p.values-c()
 for(i in 1:length(Significant[,1])){
 p.values[i]-try(t.test(positive[i,],negative[i,],alternative
 =greater)$p.value)
 }

 which(p.values0.01)


 where Significant is my matrix of  genes  and their expression in tumors
 and positive, negative are subsets of thes matrix.
 Whn p0.01, I reject the null hypothesis and I accept the alternative one,
 that I have greater gene expression in positive than in negative.
 I assume I must be doing sth wrong because the heatmap that I get with the
 genes that pass the filter of p-value is wrong.

 Could anyone help me with this?

 thanks a lot,
 Eleni



[[alternative HTML version deleted]]

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


[R] matrix inversion using solve() and matrices containing large/small values

2008-03-05 Thread gerardus vanneste
Hello

I've stumbled upon a problem for inversion of a matrix with large values,
and I haven't found a solution yet... I wondered if someone could give a
hand. (It is about automatic optimisation of a calibration process, which
involves the inverse of the information matrix)

code:

*
 macht=0.8698965
 coeff=1.106836*10^(-8)

 echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6
,-1.155094e-8,6.357603e-12)/1000
 dosis=c(0,29,70,128,201,290,396)


 dfdb -
array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8))

 dfdbtrans = aperm(dfdb)
 sigerr=sqrt(coeff*dosis^macht)
 sigmadosis = c(1:7)
 for(i in 1:7){
sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i])

}
 omega = diag(sigmadosis)
 infomatrix = dfdbtrans%*%omega%*%dfdb
**

I need the inverse of this information matrix, and

 infomatrix_inv = solve(infomatrix, tol = 10^(-43))

does not deliver adequate results (matrixproduct of infomatrix_inv and
infomatrix is not 1). Regular use of solve() delivers the error 'system is
computationally singular: reciprocal condition number: 2.949.10^(-41)'


So I went over to an eigendecomposition using eigen(): (so that infomatrix =
V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) )
in the hope this would deliver better results.)

***
 infomatrix_eigen = eigen(infomatrix)
 infomatrix_eigen_D = diag(infomatrix_eigen$values)
 infomatrix_eigen_V = infomatrix_eigen$vectors
 infomatrix_eigen_V_inv = solve(infomatrix_eigen_V)
***

however, the matrix product of these are not the same as the infomatrix
itself, only in certain parts:

 infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv
 infomatrix


Therefore, I reckon the inverse of eigendecomposition won't be correct
either.

As far as I understand, the problem is due to the very large range of data,
and therefore results in numerical problems, but I can't come up with a way
to do it otherwise.


Would anyone know how I could solve this problem?



(PS, i'm running under linux suse 10.0, latest R version with MASS libraries
(RV package))

F. Crop
UGent -- Medical Physics

[[alternative HTML version deleted]]

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


Re: [R] main title x title and y title with ggplot2

2008-03-05 Thread h . wickham
 I had the same problem and found a solution in some forums. Try this:

 p-ggplot(data, aes(x,y,fill)) + geom_point() + scale_x_continuous(your
 xlabel) + scale_y_continuous(your ylabel)

The new (more ggplot-like way) is to do:

ggplot(data, aes(x,y,fill)) + ... + opts(title = my title)

Hadley

--
http://had.co.nz/

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


Re: [R] matrix inversion using solve() and matrices containing large/small values

2008-03-05 Thread Duncan Murdoch
On 3/5/2008 8:21 AM, gerardus vanneste wrote:
 Hello
 
 I've stumbled upon a problem for inversion of a matrix with large values,
 and I haven't found a solution yet... I wondered if someone could give a
 hand. (It is about automatic optimisation of a calibration process, which
 involves the inverse of the information matrix)

If the matrix actually isn't singular, then a rescaling of the 
parameters should help a lot.  I see the diagonal of infomatrix as

  diag(infomatrix)
[1] 5.930720e-03 3.872339e+02 4.562529e+07 6.281634e+12 9.228140e+17
[6] 1.398687e+23 2.154124e+28 3.345598e+33

so multiplying the parameters by numbers on the order of the square 
roots of these entries (e.g. 10^c(-1, 1, 4, 6, 9, 12, 14, 17)), and 
redoing the rest of the calculations on that scale, should work.

Duncan Murdoch
 
 code:
 
 *
 macht=0.8698965
 coeff=1.106836*10^(-8)
 
 echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6
 ,-1.155094e-8,6.357603e-12)/1000
 dosis=c(0,29,70,128,201,290,396)
 
 
 dfdb -
 array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8))
 
 dfdbtrans = aperm(dfdb)
 sigerr=sqrt(coeff*dosis^macht)
 sigmadosis = c(1:7)
 for(i in 1:7){
 sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i])
 
 }
 omega = diag(sigmadosis)
 infomatrix = dfdbtrans%*%omega%*%dfdb
 **
 
 I need the inverse of this information matrix, and
 
 infomatrix_inv = solve(infomatrix, tol = 10^(-43))
 
 does not deliver adequate results (matrixproduct of infomatrix_inv and
 infomatrix is not 1). Regular use of solve() delivers the error 'system is
 computationally singular: reciprocal condition number: 2.949.10^(-41)'
 
 
 So I went over to an eigendecomposition using eigen(): (so that infomatrix =
 V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) )
 in the hope this would deliver better results.)
 
 ***
 infomatrix_eigen = eigen(infomatrix)
 infomatrix_eigen_D = diag(infomatrix_eigen$values)
 infomatrix_eigen_V = infomatrix_eigen$vectors
 infomatrix_eigen_V_inv = solve(infomatrix_eigen_V)
 ***
 
 however, the matrix product of these are not the same as the infomatrix
 itself, only in certain parts:
 
 infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv
 infomatrix
 
 
 Therefore, I reckon the inverse of eigendecomposition won't be correct
 either.
 
 As far as I understand, the problem is due to the very large range of data,
 and therefore results in numerical problems, but I can't come up with a way
 to do it otherwise.
 
 
 Would anyone know how I could solve this problem?
 
 
 
 (PS, i'm running under linux suse 10.0, latest R version with MASS libraries
 (RV package))
 
 F. Crop
 UGent -- Medical Physics
 
   [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


Re: [R] matrix inversion using solve() and matrices containing large/small values

2008-03-05 Thread Charilaos Skiadas
Sorry, I meant to send this to the whole list.

On Mar 5, 2008, at 8:46 AM, Charilaos Skiadas wrote:

 The problem doesn't necessarily have to do with the range of data.  
 At first level, it has to do with the simple fact that dfdb has  
 rank 6 at most, (7 at most in general, though in your case since  
 290=10*29, it is at most 6). Since matrix rank only goes down when  
 multiplying, your infomatrix is an 8x8 matrix, with rank at most 6  
 (7 if you were more lucky with that 290, still not good enough), so  
 it is certainly not invertible, even forgetting the computational  
 issues of computing the inverse.

 You would need either power smaller than 7, or a longer dosis  
 vector, I think.

 If you manage to make your problem in a case where dfdb is square,  
 then you just have to invert that, which might be easier, seeing  
 how it's a Vandermonde matrix.

 Haris Skiadas
 Department of Mathematics and Computer Science
 Hanover College

 On Mar 5, 2008, at 8:21 AM, gerardus vanneste wrote:

 Hello

 I've stumbled upon a problem for inversion of a matrix with large  
 values,
 and I haven't found a solution yet... I wondered if someone could  
 give a
 hand. (It is about automatic optimisation of a calibration  
 process, which
 involves the inverse of the information matrix)

 code:

 *
 macht=0.8698965
 coeff=1.106836*10^(-8)

 echtecoeff=c 
 (481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6
 ,-1.155094e-8,6.357603e-12)/1000
 dosis=c(0,29,70,128,201,290,396)


 dfdb -
 array(c 
 (1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7) 
 ,dim=c(7,8))

 dfdbtrans = aperm(dfdb)
 sigerr=sqrt(coeff*dosis^macht)
 sigmadosis = c(1:7)
 for(i in 1:7){
 sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^ 
 (-4),sigerr[i])

 }
 omega = diag(sigmadosis)
 infomatrix = dfdbtrans%*%omega%*%dfdb
 **

 I need the inverse of this information matrix, and

 infomatrix_inv = solve(infomatrix, tol = 10^(-43))

 does not deliver adequate results (matrixproduct of infomatrix_inv  
 and
 infomatrix is not 1). Regular use of solve() delivers the error  
 'system is
 computationally singular: reciprocal condition number: 2.949.10^ 
 (-41)'


 So I went over to an eigendecomposition using eigen(): (so that  
 infomatrix =
 V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) )
 in the hope this would deliver better results.)

 ***
 infomatrix_eigen = eigen(infomatrix)
 infomatrix_eigen_D = diag(infomatrix_eigen$values)
 infomatrix_eigen_V = infomatrix_eigen$vectors
 infomatrix_eigen_V_inv = solve(infomatrix_eigen_V)
 ***

 however, the matrix product of these are not the same as the  
 infomatrix
 itself, only in certain parts:

 infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv
 infomatrix


 Therefore, I reckon the inverse of eigendecomposition won't be  
 correct
 either.

 As far as I understand, the problem is due to the very large range  
 of data,
 and therefore results in numerical problems, but I can't come up  
 with a way
 to do it otherwise.


 Would anyone know how I could solve this problem?



 (PS, i'm running under linux suse 10.0, latest R version with MASS  
 libraries
 (RV package))

 F. Crop
 UGent -- Medical Physics

  [[alternative HTML version deleted]]






Haris Skiadas
Department of Mathematics and Computer Science
Hanover College

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


Re: [R] t.test p-Value

2008-03-05 Thread Eleni Christodoulou
On Wed, Mar 5, 2008 at 2:05 PM, ian white [EMAIL PROTECTED] wrote:

 Don't you need to make some allowance for multiple testing? E.g. to get
 a experiment-wise significance level of 0.01 you need

 which(p.values  very small number)

 where the very small number is approximately 0.01/(total number of
 genes).

 On Wed, 2008-03-05 at 14:38 +0200, Eleni Christodoulou wrote:
  I am sorry, the test is unpaired...But my question remains
 
  Thanks,
  Eleni
 
  On Wed, Mar 5, 2008 at 2:33 PM, Eleni Christodoulou [EMAIL PROTECTED]
 
  wrote:
 
   Hello list,
  
   I am trying to apply the paired t.test between diseased and not
 diseased
   patients to identify genes that are more expressed in the one
 situation
   under the other. In order to retrieve the genes that are more
 expressed in
   the positive disease state I do:
   p.values-c()
   for(i in 1:length(Significant[,1])){
   p.values[i]-try(t.test(positive[i,],negative[i,],alternative
   =greater)$p.value)
   }
  
   which(p.values0.01)
  
  
   where Significant is my matrix of  genes  and their expression in
 tumors
   and positive, negative are subsets of thes matrix.
   Whn p0.01, I reject the null hypothesis and I accept the alternative
 one,
   that I have greater gene expression in positive than in negative.
   I assume I must be doing sth wrong because the heatmap that I get with
 the
   genes that pass the filter of p-value is wrong.
  
   Could anyone help me with this?
  
   thanks a lot,
   Eleni
  
  
 
[[alternative HTML version deleted]]
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 



[[alternative HTML version deleted]]

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


[R] nls: different results if applied to normal or linearized data

2008-03-05 Thread Wolfgang Waser
Dear all,

I did a non-linear least square model fit

y ~ a * x^b 

(a)  nls(y ~ a * x^b, start=list(a=1,b=1))

to obtain the coefficients a  b.

I did the same with the linearized formula, including a linear model

log(y) ~ log(a) + b * log(x)

(b)  nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1))
(c)  lm(log10(y) ~ log10(x))

I expected coefficient b to be identical for all three cases. Hoever, using my 
dataset, coefficient b was:
(a) 0.912
(b) 0.9794
(c) 0.9794

Coefficient a also varied between option (a) and (b), 107.2 and 94.7, 
respectively.

Is this supposed to happen? Which is the correct coefficient b?


Regards,

Wolfgang

-- 
Laboratory of Animal Physiology
Department of Biology
University of Turku
FIN-20014 Turku
Finland

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


Re: [R] Asking, are simple effects different from 0

2008-03-05 Thread Chuck Cleland
On 3/4/2008 2:45 PM, Jarrett Byrnes wrote:
 Hello, R-i-zens.  I'm working on an data set with a factorial ANOVA  
 that has a significant interaction.  I'm interested in seeing whether  
 the simple effects are different from 0, and I'm pondering how to do  
 this.  So, I have
 
 my.anova-lm(response ~ trtA*trtB)
 
 The output for which gives me a number of coefficients and whether  
 they are different from 0.  However, I want the simple effects only,  
 incorporating their intercepts, with their error estimates.  Can I  
 somehow manipulate this object to get that?  Or, would a good shortcut  
 be
 
 my.simple.anova-lm(response ~ trtA:trtB + 0)
 
 and then use those coefficients and their error estimates?

   One approach would be to use glht() in the multcomp package.  You 
need to work out how to formulate the matrix of coefficients that give 
the desired contrasts.  Here is an example using the warpbreaks data frame:

fm - lm(breaks ~ tension*wool, data=warpbreaks)

# names(coef(fm))
# (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB

cm - rbind(
A vs. B at L = c(0, 0, 0,-1, 0, 0),
A vs. B at M = c(0, 0, 0,-1,-1, 0),
A vs. B at H = c(0, 0, 0,-1, 0,-1),
M vs. L at A = c(0, 1, 0, 0, 0, 0),
M vs. H at A = c(0, 1,-1, 0, 0, 0),
L vs. H at A = c(0, 0,-1, 0, 0, 0),
M vs. L at B = c(0, 1, 0, 0, 1, 0),
M vs. H at B = c(0, 1,-1, 0, 1,-1),
L vs. H at B = c(0, 0,-1, 0, 0,-1))

library(multcomp)

summary(glht(fm, linfct = cm), test = adjusted(type=none))

  Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks)

Linear Hypotheses:
   Estimate Std. Error t value  p value
A vs. B at L == 0  16. 5.1573   3.167 0.002677 **
A vs. B at M == 0  -4.7778 5.1573  -0.926 0.358867
A vs. B at H == 0   5.7778 5.1573   1.120 0.268156
M vs. L at A == 0 -20.5556 5.1573  -3.986 0.000228 ***
M vs. H at A == 0  -0.5556 5.1573  -0.108 0.914665
L vs. H at A == 0  20. 5.1573   3.878 0.000320 ***
M vs. L at B == 0   0.5556 5.1573   0.108 0.914665
M vs. H at B == 0  10. 5.1573   1.939 0.058392 .
L vs. H at B == 0   9. 5.1573   1.831 0.073270 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- none method)

 If so, I realize that R gives me t tests for each coefficient.  Now,  
 for those I know I'm using the residual degrees of freedom.  Would it  
 then be more appropriate to use those, all with the same residual DF  
 and apply a bonferroni correction, or, use the mean and SE estimate  
 with the sample size for that particular treatment and perform an  
 uncorrected one sample t-test to see if the value is different from 0?

   I won't comment on whether to adjust and if so how, but you can 
implement various adjustments when summarizing.  For example:

summary(glht(fm, linfct = cm), test = adjusted(type=bonferroni))

  Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks)

Linear Hypotheses:
   Estimate Std. Error t value p value
A vs. B at L == 0  16. 5.1573   3.167 0.02409 *
A vs. B at M == 0  -4.7778 5.1573  -0.926 1.0
A vs. B at H == 0   5.7778 5.1573   1.120 1.0
M vs. L at A == 0 -20.5556 5.1573  -3.986 0.00205 **
M vs. H at A == 0  -0.5556 5.1573  -0.108 1.0
L vs. H at A == 0  20. 5.1573   3.878 0.00288 **
M vs. L at B == 0   0.5556 5.1573   0.108 1.0
M vs. H at B == 0  10. 5.1573   1.939 0.52553
L vs. H at B == 0   9. 5.1573   1.831 0.65943
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- bonferroni method)

 Sorry for the bonehead question, but it's a situation I haven't seen  
 before.
 
 -Jarrett
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code. 

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


[R] (no subject)

2008-03-05 Thread Zhao, Wei (Cancer Center)
Dear all,

 

I try to import a SPSS.por dataset with about 6000 cases and 650
variables with R commander 

and got error messages:

'Error: /temp/t.por  use.value.labels=true,
max.value.label=inf,

'Error: data is not data frame and cannot be attached'

 

Any comment?

 

Thanks in advance.

 

wz

http://geocities.com/wendyzhao78/HW.html 

 

 


[[alternative HTML version deleted]]

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


Re: [R] nls: different results if applied to normal or linearized data

2008-03-05 Thread Gabor Grothendieck
Write out the objective functions that they are minimizing and it
will be clear they are different so you can't expect the same
results.

On Wed, Mar 5, 2008 at 8:53 AM, Wolfgang Waser [EMAIL PROTECTED] wrote:
 Dear all,

 I did a non-linear least square model fit

 y ~ a * x^b

 (a)  nls(y ~ a * x^b, start=list(a=1,b=1))

 to obtain the coefficients a  b.

 I did the same with the linearized formula, including a linear model

 log(y) ~ log(a) + b * log(x)

 (b)  nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1))
 (c)  lm(log10(y) ~ log10(x))

 I expected coefficient b to be identical for all three cases. Hoever, using my
 dataset, coefficient b was:
 (a) 0.912
 (b) 0.9794
 (c) 0.9794

 Coefficient a also varied between option (a) and (b), 107.2 and 94.7,
 respectively.

 Is this supposed to happen? Which is the correct coefficient b?


 Regards,

 Wolfgang

 --
 Laboratory of Animal Physiology
 Department of Biology
 University of Turku
 FIN-20014 Turku
 Finland

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


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


[R] Updating packages from one hard-drive to another, after upgrade of R

2008-03-05 Thread Bob Green

The CPU on my computer 'died' and I have had to purchase a new 
computer. I have just installed v 2.6.2. My previous computer had v 
2.5.1 and a large number of files in the library folder. These files 
have been copied to a partition on my new hard drive, along with the 
old R installation.

Can I just copy all the folders/files in the old library to the new 
v2.6.2 library or do the packages need to be downloaded again? There 
are quite a number of folders so it would be time-consuming to have 
to do a fresh download.

Any assistance is appreciated,

Bob

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


[R] Import SPSS.por file

2008-03-05 Thread Zhao, Wei (Cancer Center)
Dear all,

 

I try to import a SPSS.por dataset with about 6000 cases and 650
variables with R commander 

and got error messages:

'Error: /temp/t.por  use.value.labels=true,
max.value.label=inf,

'Error: data is not data frame and cannot be attached'

 

Any comment?

 

Thanks in advance.

 

wz

http://geocities.com/wendyzhao78/HW.html 

 

 


[[alternative HTML version deleted]]

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


Re: [R] coxme - fitting random treatment effect nested within centre

2008-03-05 Thread Terry Therneau
 included message 
 Dear all,

I am using coxme function in Kinship library to fit random treatment effect 
nested within centre. I got 3 treatments (0,1,2) and 3 centres. I used 
following 
commands, but got an error. 

 ugroup=paste(rep(1:3,each=3),rep(0:2,3),sep='/')
 
mat1=bdsmatrix(rep(c(1,1,1,1,1,1,1,1,1),3),blocksize=rep(3,3),dimnames=list(ugro
up,ugroup))
 
mat2=bdsmatrix(rep(c(0,0,0,0,0,0,0,0,1),3),blocksize=rep(3,3),dimnames=list(ugro
up,ugroup))
 group=paste(dat1$centre,dat1$treat,sep='/')
 coxme(Surv(time,status) ~ as.factor(treat), data=dat1,random= 
~1|group,varlist=list(mat1,mat2),rescale=F,pdcheck=FALSE)

Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, : 
Random effects variance is not spd
 
Could anyone help me correcting this error? 
Many thanks in advance.
Ruwanthi

--- end inclusion -

 The build your own bdsmatrix style in coxme is for variance structures that 
are not built in, like pedigree data.  Your problem is one that coxme can do 
directly, so it is easier to call the routine simply:
 
  fit - coxme(Surv(time, status) ~ factor(treat), data=data1,
 random = ~1 | centre/treat)
   
Terry Therneau

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


[R] legend for several graphics

2008-03-05 Thread Georg Otto
Hi,

I am trying to generate a figure of 9 plots that are contained in one
device by using

par(mfrow = c(3,3,))

I would like to have 1 common legend for all 9 plots somewhere outside
of the plotting area (as opposed to one legend inside each of the 9
plots, which the function legend() seems to generate by default).

Any hint how to do this?

Best,

Georg

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


Re: [R] nls: different results if applied to normal or linearized data

2008-03-05 Thread Martin Elff
On Wednesday 05 March 2008 (14:53:27), Wolfgang Waser wrote:
 Dear all,

 I did a non-linear least square model fit

 y ~ a * x^b

 (a)  nls(y ~ a * x^b, start=list(a=1,b=1))

 to obtain the coefficients a  b.

 I did the same with the linearized formula, including a linear model

 log(y) ~ log(a) + b * log(x)

 (b)  nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1))
 (c)  lm(log10(y) ~ log10(x))

 I expected coefficient b to be identical for all three cases. Hoever, using
 my dataset, coefficient b was:
 (a) 0.912
 (b) 0.9794
 (c) 0.9794

 Coefficient a also varied between option (a) and (b), 107.2 and 94.7,
 respectively.

Models (a) and (b) entail different distributions of the dependent variable y
and different ranges of values that y may take.
(a) implies that y has, conditionally on x, a normal distribution and
has a range of feasible values from -Inf to +Inf.
(b) and (c) imply that log(y) has a normal distribution, that is, 
y has a log-normal distribution and can take values from zero to +Inf.

 Is this supposed to happen? 
Given the above considerations, different results with respect to the 
intercept are definitely to be expected.

 Which is the correct coefficient b? 
That depends - is y strictly non-negative or not ...

Just my 20 cents...

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


Re: [R] Updating packages from one hard-drive to another, after upgrade of R

2008-03-05 Thread ONKELINX, Thierry
Bob,

You can copy the files from the packages to your new computer. Then run
update.packages(checkBuilt = TRUE).

That should do.

HTH,

Thierry




ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium 
tel. + 32 54/436 185
[EMAIL PROTECTED] 
www.inbo.be 

Do not put your faith in what statistics say until you have carefully
considered what they do not say.  ~William W. Watt
A statistical analysis, properly conducted, is a delicate dissection of
uncertainties, a surgery of suppositions. ~M.J.Moroney

-Oorspronkelijk bericht-
Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
Namens Bob Green
Verzonden: woensdag 5 maart 2008 15:07
Aan: r-help@r-project.org
Onderwerp: [R] Updating packages from one hard-drive to another,after
upgrade of R


The CPU on my computer 'died' and I have had to purchase a new 
computer. I have just installed v 2.6.2. My previous computer had v 
2.5.1 and a large number of files in the library folder. These files 
have been copied to a partition on my new hard drive, along with the 
old R installation.

Can I just copy all the folders/files in the old library to the new 
v2.6.2 library or do the packages need to be downloaded again? There 
are quite a number of folders so it would be time-consuming to have 
to do a fresh download.

Any assistance is appreciated,

Bob

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

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


Re: [R] Updating packages from one hard-drive to another, after upgrade of R

2008-03-05 Thread Gabor Grothendieck
See ?.libPaths for info on setting the location of your libraries.
If you do want to copy them the copydir.bat utility in
 batchfiles.googlecode.com
can do that.

On Wed, Mar 5, 2008 at 9:07 AM, Bob Green [EMAIL PROTECTED] wrote:

 The CPU on my computer 'died' and I have had to purchase a new
 computer. I have just installed v 2.6.2. My previous computer had v
 2.5.1 and a large number of files in the library folder. These files
 have been copied to a partition on my new hard drive, along with the
 old R installation.

 Can I just copy all the folders/files in the old library to the new
 v2.6.2 library or do the packages need to be downloaded again? There
 are quite a number of folders so it would be time-consuming to have
 to do a fresh download.

 Any assistance is appreciated,

 Bob

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


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


Re: [R] matrix inversion using solve() and matrices containing large/small values

2008-03-05 Thread Douglas Bates
On Wed, Mar 5, 2008 at 7:43 AM, Duncan Murdoch [EMAIL PROTECTED] wrote:
 On 3/5/2008 8:21 AM, gerardus vanneste wrote:
   Hello
  
   I've stumbled upon a problem for inversion of a matrix with large values,
   and I haven't found a solution yet...

Someone with experience in numerical linear algebra will immediately
focus on the words inversion of a matrix in that statement.

There is a sort of a meta-result in numerical linear algebra that if
the best way you can formulate an answer to your problem involves
inversion of a matrix (without specifying a special structure like
diagonal or bidiagonal or unit triangular or ...) then you need to
think about your problem in more detail.

Certainly a formula may be written out in terms of the inverse of an
information matrix but it is a bad idea to try to calculate the result
that way.  If you have an information matrix is it positive definite?
If so, you should be working with the Cholesky decomposition of the
matrix or perhaps the QR decomposition of a model matrix that
generates the information matrix.

Think in terms of the factors of a matrix and how you can reexpress
the calculation using them.

   I wondered if someone could give a
   hand. (It is about automatic optimisation of a calibration process, which
   involves the inverse of the information matrix)

  If the matrix actually isn't singular, then a rescaling of the
  parameters should help a lot.  I see the diagonal of infomatrix as

diag(infomatrix)
  [1] 5.930720e-03 3.872339e+02 4.562529e+07 6.281634e+12 9.228140e+17
  [6] 1.398687e+23 2.154124e+28 3.345598e+33

  so multiplying the parameters by numbers on the order of the square
  roots of these entries (e.g. 10^c(-1, 1, 4, 6, 9, 12, 14, 17)), and
  redoing the rest of the calculations on that scale, should work.

  Duncan Murdoch


 
   code:
  
   *
   macht=0.8698965
   coeff=1.106836*10^(-8)
  
   echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6
   ,-1.155094e-8,6.357603e-12)/1000
   dosis=c(0,29,70,128,201,290,396)
  
  
   dfdb -
   
 array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8))
  
   dfdbtrans = aperm(dfdb)
   sigerr=sqrt(coeff*dosis^macht)
   sigmadosis = c(1:7)
   for(i in 1:7){
   
 sigmadosis[i]=ifelse(sigerr[i]2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i])
  
   }
   omega = diag(sigmadosis)
   infomatrix = dfdbtrans%*%omega%*%dfdb
   **
  
   I need the inverse of this information matrix, and
  
   infomatrix_inv = solve(infomatrix, tol = 10^(-43))
  
   does not deliver adequate results (matrixproduct of infomatrix_inv and
   infomatrix is not 1). Regular use of solve() delivers the error 'system is
   computationally singular: reciprocal condition number: 2.949.10^(-41)'
  
  
   So I went over to an eigendecomposition using eigen(): (so that infomatrix 
 =
   V D V^(-1) == infomatrix^(-1)= V D^(-1) V^(-1) )
   in the hope this would deliver better results.)
  
   ***
   infomatrix_eigen = eigen(infomatrix)
   infomatrix_eigen_D = diag(infomatrix_eigen$values)
   infomatrix_eigen_V = infomatrix_eigen$vectors
   infomatrix_eigen_V_inv = solve(infomatrix_eigen_V)
   ***
  
   however, the matrix product of these are not the same as the infomatrix
   itself, only in certain parts:
  
   infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv
   infomatrix
  
  
   Therefore, I reckon the inverse of eigendecomposition won't be correct
   either.
  
   As far as I understand, the problem is due to the very large range of data,
   and therefore results in numerical problems, but I can't come up with a way
   to do it otherwise.
  
  
   Would anyone know how I could solve this problem?
  
  
  
   (PS, i'm running under linux suse 10.0, latest R version with MASS 
 libraries
   (RV package))
  
   F. Crop
   UGent -- Medical Physics
  
 [[alternative HTML version deleted]]
  
   __
   R-help@r-project.org mailing list
   https://stat.ethz.ch/mailman/listinfo/r-help
   PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
   and provide commented, minimal, self-contained, reproducible code.

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


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


[R] regex sulotion for seperating number and string

2008-03-05 Thread sun
I have strings contain postcode and letters, some seperated with blank, some 
with comma, and some hasn't seperated. eg, 2324gz  2567 HK 3741,BF

I want to seperate the number and letters into two new variables.

I know this should be quite basic question, but searched on regex syntax and 
that seems a bit scarey to me, any one can shot me a quick solution on this 
particular question?

thanks,
Sun

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


Re: [R] vertex labels in igraph from adjacency matrix

2008-03-05 Thread Mark W Kimpel
Looks like I turned an off my one error into an off by two error by 
adding rather than subtracting. Clearly a logic error on my part.

Also, which.max is clearly superior as it results in half as many 
function calls.

Thanks guys!

As an aside, although igraph may use the C indexing convention, R users 
and their code, is much more in tune with indexing beginning at 1. 
Adjusting this in the R interface to igraph may solve many headaches 
down the road.

Mark

Mark W. Kimpel MD  ** Neuroinformatics ** Dept. of Psychiatry
Indiana University School of Medicine

15032 Hunter Court, Westfield, IN  46074

(317) 490-5129 Work,  Mobile  VoiceMail
(317) 204-4202 Home (no voice mail please)

mwkimpelatgmaildotcom

**


Gabor Csardi wrote:
 On Wed, Mar 05, 2008 at 02:27:21AM -0500, Charilaos Skiadas wrote:
 [...]
 Btw, you will likely want to take the betweenness call out, and call  
 it once and store the result, instead of calling it twice (well,  
 assuming the graph is largish). Or even better, use which.max:

 which.max(betweenness(graph = my.graph, v=V(my.graph), directed =  
 FALSE))
 
 This is almost good, but there is a catch, in igraph vertices are 
 numbered from zero. So if you want an igraph vertex id, then you 
 need to subtract one from this, i.e.:
 
 maxb - which.max(betweennness(my.graph, directed=FALSE))-1
 
 You can double check it:
 
 betweenness(my.graph, maxb, directed=FALSE)
 
 Gabor
 
 PS. there is also an igraph mailing list, see the igraph homepage
 at igraph.sf.net
 
 Haris Skiadas
 Department of Mathematics and Computer Science
 Hanover College

 [...]


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


[R] problem with geepack

2008-03-05 Thread Giacomo Santini
Hi all

I am analyzing a data set containing information about the behaviour of 
marine molluscs on a vertical wall. Since I have replicate observations 
on the same individuals I was thinking to use the geepack library.

The data are organised in a dataframe with the following variables

Date = date of sampling,
Size = dimensions (mm)
Activity duration of activity (min)
Water = duration of splashing by waves
Hgt = resting eight of each specimen before activity begin
Individual = a code indicating the id of the specimen.

I have up to 12 replicate observations for individual. Some observation 
are missing and I organized the data frame to have exactly 12 rows for 
each specimen, with NAs where there is a missing observation.

The following model worked fine:

gee1-geese(Activity~Water, id=Individual, data=dataF, family=gaussian)

but when I use other variables e.g

gee2-geese(Activity~Hgt+Size+Water, id=Individual, data=dataF, 
family=gaussian)

I get the error message

Error in geese.fit(x, y, id, offset, soffset, w, waves, zsca, zcor, 
corp,  :
  nrow(zsca) and length(y) not match

which I am not able to understand.

The same problem has been reported in the list in 2006, but I have not 
found any response to it.

Any suggestion?


Giacomo

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


[R] Constrained regression

2008-03-05 Thread Carlos Alzola
I would like to acknowledge the answers I received from Tom Filloon, Mike
Cheung and Berwyn Turlach.
Berwyn's response was exactly what I needed. Use solve.QP from the quadprog
package in R. S-Plus has the equivalent function solveQP in the NuOpt
module.

Berwyn's response is below

G'day Carlos,

On Mon, Mar 3, 2008 at 11:52 AM
Carlos Alzola [EMAIL PROTECTED] wrote:

  I am trying to get information on how to fit a linear regression with 
 constrained parameters. Specifically, I have 8 predictors , their 
 coeffiecients should all be non-negative and add up to 1. I understand 
 it is a quadratic programming problem but I have no experience in the 
 subject. I searched the archives but the results were inconclusive.

  Could someone provide suggestions and references to the literature, 
 please?

A suggestion:

 library(MASS)   ## to access the Boston data
 designmat - model.matrix(medv~., data=Boston) Dmat - 
 crossprod(designmat, designmat) dvec - crossprod(designmat, 
 Boston$medv) Amat - cbind(1, diag(NROW(Dmat))) bvec - c(1, 
 rep(0,NROW(Dmat)) meq - 1
 library(quadprog)
 res - solve.QP(Dmat, dvec, Amat, bvec, meq)

The solution seems to contain values that are, for all practical purposes,
actually zero:

 res$solution
 [1]  4.535581e-16  2.661931e-18  1.016929e-01 -1.850699e-17  [5]
1.458219e-16 -3.892418e-15  8.544939e-01  0.00e+00  [9]  2.410742e-16
2.905722e-17 -5.700600e-20 -4.227261e-17 [13]  4.381328e-02 -3.723065e-18

So perhaps better:

 zapsmall(res$solution)
 [1] 0.000 0.000 0.1016929 0.000 0.000 0.000  [7]
0.8544939 0.000 0.000 0.000 0.000 0.000 [13] 0.0438133
0.000

So the estimates seem to follow the constraints.

And the unconstrained solution is:

 res$unconstrainted.solution
 [1]  3.645949e+01 -1.080114e-01  4.642046e-02  2.055863e-02  [5]
2.686734e+00 -1.776661e+01  3.809865e+00  6.922246e-04  [9] -1.475567e+00
3.060495e-01 -1.233459e-02 -9.527472e-01 [13]  9.311683e-03 -5.247584e-01

which seems to coincide with what lm() thinks it should be:

 coef(lm(medv~., Boston))
  (Intercept)  crimzn indus  chas 
 3.645949e+01 -1.080114e-01  4.642046e-02  2.055863e-02  2.686734e+00 
  noxrm   age   dis   rad 
-1.776661e+01  3.809865e+00  6.922246e-04 -1.475567e+00  3.060495e-01 
  tax   ptratio black lstat 
-1.233459e-02 -9.527472e-01  9.311683e-03 -5.247584e-01 

So there seem to be no numeric problems.  Otherwise we could have done
something else (e.g calculate the QR factorization of the design matrix, say
X, and give the R factor to solve.QP, instead of calculating X'X and giving
that one to solve.QP).

If the intercept is not supposed to be included in the set of constrained
estimates, then something like the following can be done:

 Amat[1,] - 0
 res - solve.QP(Dmat, dvec, Amat, bvec, meq)
 zapsmall(res$solution)
 [1] 6.073972 0.00 0.109124 0.00 0.00 0.00 0.863421  [8]
0.00 0.00 0.00 0.00 0.00 0.027455 0.00

Of course, since after the first command in that last block the second
column of Amat contains only zeros
 Amat[,2]
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
we might as well have removed it (and the corresponding entry in bvec)
 Amat - Amat[, -2]
 bvec - bvec[-2]
before calling solve.QP().

Note, the Boston data set was only used to illustrate how to fit such
models, I do not want to imply that these models are sensible for these
data. :-)

Hope this helps.

Cheers,

Berwin


Carlos Alzola
[EMAIL PROTECTED]
(703) 242-6747 


[[alternative HTML version deleted]]

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


[R] XLSolutions 9 Courses: Upcoming March-April 2008 R/S+ Course Schedule by XLSolutions Corp

2008-03-05 Thread [EMAIL PROTECTED]
Our March-April 2008 R/S+ course schedule is now available. Please check
out this link for additional information  and direct enquiries to Sue
Turner [EMAIL PROTECTED]  Phone: 206 686 1578
 
 --Can't see your city?   Please email us!  --
      Ask for  Group Discount ---

www.xlsolutions-corp.com/courselist.htm
 
(1) R/S System: Advanced Programming 
 
*** San Francisco / April 28-29, 2008 *** 
*** Salt Lake City / April 21-22, 2008 *** 

(2) R/S Fundamentals and Programming Techniques 
 
*** San Francisco / March 27-28, 2008  ***
*** New York City / March 25-26, 2008  ***

*** Seattle /  April 21-22
 
(3) Data Mining: Practical Tools and Techniques in R/Splus  
 
*** San Franciso / April 24-25, 2008 ***
 

(4) R/S System: Complementing and Extending Statistical Computing for
SAS Users 
*** Raleigh / April 28-29, 2008 ***  
 
(5) Introduction to R/S+ programming: Microarrays Analysis and
Bioconductor.
*** Los Angeles   /  April 25-26, 2008*** 
 
(6) Microarrays Data Analysis with R/S+ and GGobi  
*** New York City /  April 28-29, 2008*** 
 
 
www.xlsolutions-corp.com/courselist.htm
 
 

 Payment due AFTER the class
 Email us for group discounts. 
 Email Sue Turner: [EMAIL PROTECTED] 
 Phone: 206-686-1578 
 Visit us: www.xlsolutions-corp.com/courselist.htm 
 Please let us know if you and your colleagues are interested in this 
 class to take advantage of group discount. Register now to secure your 
 seat! 
 
 Cheers, 
 Elvis Miller, PhD 
 Manager Training. 
 XLSolutions Corporation 
 206 686 1578 
 www.xlsolutions-corp.com

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


Re: [R] regex sulotion for seperating number and string

2008-03-05 Thread vincenzo . 2 . di-iorio
Hi Sun,

vec -  c(2324gz,2567 HK,3741,BF)

vec1 - gsub('[^[:digit:]]','',vec)

vec2 - gsub('[^[:alpha:]]','',vec)
 
 vec1
[1] 2324 2567 3741
 vec2
[1] gz HK BF


Cheers
Vincenzo


---
Vincenzo Luca Di Iorio
Consultant PME User support - GSK RD Limited
---




sun [EMAIL PROTECTED] 
Sent by: [EMAIL PROTECTED]
05-Mar-2008 15:51
 
To
[EMAIL PROTECTED]
cc

Subject
[R] regex sulotion for seperating number and string






I have strings contain postcode and letters, some seperated with blank, 
some 
with comma, and some hasn't seperated. eg, 2324gz  2567 HK 3741,BF

I want to seperate the number and letters into two new variables.

I know this should be quite basic question, but searched on regex syntax 
and 
that seems a bit scarey to me, any one can shot me a quick solution on 
this 
particular question?

thanks,
Sun

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



[[alternative HTML version deleted]]

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


[R] CROSSOVER TRIALS IN R (Binary Outcomes)

2008-03-05 Thread Boikanyo Makubate
I will like to analyse a binary cross over design using the random 
effects model. The probability of success is assumed to be logistic. 
Suppose as an example, we have 4 subjects undergoing a crossover design, 
where the outcome is either success or failure. The first two subjects 
receive treatment A first followed by treatment B. The remaining two 
subjects receive treatments in the reverse order. The outcomes for the 
subjects is the sequence AB are as follows: (0,1) and (0,0). While the 
outcomes for the subjects in the sequence BA are (1,1) and (1,0).

How can i analyse this using R. I have done the problem with PROC 
NLMIXED in SAS, I simply want to compare the results from SAS with those 
from R. Please help. R-Codes will be highly appreciated.

Boikanyo.

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


Re: [R] Asking, are simple effects different from 0

2008-03-05 Thread jebyrnes

Huh.  Very interesting.  I haven't really worked with manipulating contrast
matrices before, save to do a prior contrasts.  Could you explain the matrix
you laid out just a bit more so that I can generalize it to my case?  


Chuck Cleland wrote:
 
 
One approach would be to use glht() in the multcomp package.  You 
 need to work out how to formulate the matrix of coefficients that give 
 the desired contrasts.  Here is an example using the warpbreaks data
 frame:
 
 fm - lm(breaks ~ tension*wool, data=warpbreaks)
 
 # names(coef(fm))
 # (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB
 
 cm - rbind(
 A vs. B at L = c(0, 0, 0,-1, 0, 0),
 A vs. B at M = c(0, 0, 0,-1,-1, 0),
 A vs. B at H = c(0, 0, 0,-1, 0,-1),
 M vs. L at A = c(0, 1, 0, 0, 0, 0),
 M vs. H at A = c(0, 1,-1, 0, 0, 0),
 L vs. H at A = c(0, 0,-1, 0, 0, 0),
 M vs. L at B = c(0, 1, 0, 0, 1, 0),
 M vs. H at B = c(0, 1,-1, 0, 1,-1),
 L vs. H at B = c(0, 0,-1, 0, 0,-1))
 
 library(multcomp)
 
 summary(glht(fm, linfct = cm), test = adjusted(type=none))
 
   Simultaneous Tests for General Linear Hypotheses
 
 Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks)
 
 Linear Hypotheses:
Estimate Std. Error t value  p value
 A vs. B at L == 0  16. 5.1573   3.167 0.002677 **
 A vs. B at M == 0  -4.7778 5.1573  -0.926 0.358867
 A vs. B at H == 0   5.7778 5.1573   1.120 0.268156
 M vs. L at A == 0 -20.5556 5.1573  -3.986 0.000228 ***
 M vs. H at A == 0  -0.5556 5.1573  -0.108 0.914665
 L vs. H at A == 0  20. 5.1573   3.878 0.000320 ***
 M vs. L at B == 0   0.5556 5.1573   0.108 0.914665
 M vs. H at B == 0  10. 5.1573   1.939 0.058392 .
 L vs. H at B == 0   9. 5.1573   1.831 0.073270 .
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 (Adjusted p values reported -- none method)
 
 
 

-- 
View this message in context: 
http://www.nabble.com/Asking%2C-are-simple-effects-different-from-0-tp15835552p15852362.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] loop

2008-03-05 Thread Neuer Arkadasch
Hello all,
   
  I am trying to use 
   
  m - seq(-1,1,0.1)
  x1 - vector()
  x2 - vector()
  for(i in m){
  x1[i] - i
  x2[i] - i^2
  }
  dat - data.frame(x1,x2)
  But, I have  false result
  dat
x1 x2
  1 1  1
   
  could some tell me how it is possible to do this?
   
  Thank you!
   
   

   
-

[[alternative HTML version deleted]]

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


[R] rrp.impute: for what sizes does it work?

2008-03-05 Thread Werner Wernersen
Hi,

I have a survey dataset of about 2 observations
where for 2 factor variables I have about 200 missing
values each. I want to impute these using 10 possibly
explanatory variables which are a mixture of integers
and factors.

Since I was quite intrigued by the concept of rrp I
wanted to use it but it takes ages and terminates with
an error. First time it stopped complaining about too
little memory which I increased then from 1.5 to 2GB. 
Now it terminates after a long time with this:
Error in nn[[i]] : subscript out of bounds

Has anybody encountered this problem before, is it due
to my data? Could you recommend another package for
such imputation? 
I looked at the aregImpute in Hmisc but I don't really
understand what it is doing.

Thanks a lot,
 Werner


  E-Mails jetzt auf Ihrem Handy.

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


[R] Need help for calculating cross-correlation between 4 multivariate time series data

2008-03-05 Thread sandeep jose joseph
Hi all,
 I would like to know whether there is any function in R were i can
find the cross-correlation of two or more multivariate (time series) data. I
tried the function ccf() but it seems like to have two univariate datasets.
Please let me know.
sincerely,
sandeep

-- 
Sandeep Joseph PhD
Post Doctoral Associate
Center for Tropical  Emerging Global Diseases
Paul D. Coverdell Center, Rm 335A
500 D. W. Brooks Drive
Athens, GA 30602-7394
ph 404-578-6790

2091 south Miledge Ave
Apt A6, Athens,
GA-30605.
phone: 706-850-0056

[[alternative HTML version deleted]]

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


Re: [R] nls: different results if applied to normal or linearized data

2008-03-05 Thread Prof Brian Ripley
On Wed, 5 Mar 2008, Wolfgang Waser wrote:

 Dear all,

 I did a non-linear least square model fit

 y ~ a * x^b

 (a)  nls(y ~ a * x^b, start=list(a=1,b=1))

 to obtain the coefficients a  b.

 I did the same with the linearized formula, including a linear model

 log(y) ~ log(a) + b * log(x)

 (b)  nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1))
 (c)  lm(log10(y) ~ log10(x))

 I expected coefficient b to be identical for all three cases. Hoever, using my
 dataset, coefficient b was:
 (a) 0.912
 (b) 0.9794
 (c) 0.9794

 Coefficient a also varied between option (a) and (b), 107.2 and 94.7,
 respectively.

 Is this supposed to happen? Which is the correct coefficient b?

Yes.  You are fitting by least-squares on two different scales: 
differences in y and differences in log(y) are not comparable.

Both are correct solutions to different problems.  Since we have no idea 
what 'x' and 'y' are, we cannot even guess which is more appropriate in 
your context.



 Regards,

 Wolfgang

 --
 Laboratory of Animal Physiology
 Department of Biology
 University of Turku
 FIN-20014 Turku
 Finland

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


-- 
Brian D. Ripley,  [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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] loop

2008-03-05 Thread Giovanni Petris

 Date: Wed, 05 Mar 2008 15:59:59 +0100 (CET)
 From: Neuer Arkadasch [EMAIL PROTECTED]
 Sender: [EMAIL PROTECTED]
 Precedence: list
 
 Hello all,

   I am trying to use 

   m - seq(-1,1,0.1)
   x1 - vector()
   x2 - vector()
   for(i in m){
   x1[i] - i
   x2[i] - i^2
   }
   dat - data.frame(x1,x2)
   But, I have  false result
   dat
 x1 x2
   1 1  1


It does not seem to be false: I am getting the same result. You should
probably explain what you expected and what you are trying to achieve.

Giovanni

   could some tell me how it is possible to do this?

   Thank you!


 

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

-- 

Giovanni Petris  [EMAIL PROTECTED]
Associate Professor
Department of Mathematical Sciences
University of Arkansas - Fayetteville, AR 72701
Ph: (479) 575-6324, 575-8630 (fax)
http://definetti.uark.edu/~gpetris/

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


Re: [R] loop

2008-03-05 Thread Ravi Varadhan
Try this:

m - seq(-1,1,0.1)
  x1 - vector(length=length(m))
  x2 - vector(length=length(m))
  for(i in m){
  x1[i] - i
  x2[i] - i^2
  }
  dat - data.frame(x1,x2)

Ravi.


---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: [EMAIL PROTECTED]

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

 





-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Neuer Arkadasch
Sent: Wednesday, March 05, 2008 10:00 AM
To: [EMAIL PROTECTED]
Subject: [R] loop

Hello all,
   
  I am trying to use 
   
  m - seq(-1,1,0.1)
  x1 - vector()
  x2 - vector()
  for(i in m){
  x1[i] - i
  x2[i] - i^2
  }
  dat - data.frame(x1,x2)
  But, I have  false result
  dat
x1 x2
  1 1  1
   
  could some tell me how it is possible to do this?
   
  Thank you!
   
   

   
-

[[alternative HTML version deleted]]

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

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


Re: [R] loop

2008-03-05 Thread Yinghai Deng
m - seq(-1,1,0.1)
  x1 - c()
  x2 - c()
  for(i in 1:length(m)){
  x1[i] - m[i]
  x2[i] - m[i]^2
  }
  dat - data.frame(x1,x2)

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Behalf Of Neuer Arkadasch
Sent: March 5, 2008 10:00 AM
To: [EMAIL PROTECTED]
Subject: [R] loop


Hello all,

  I am trying to use

  m - seq(-1,1,0.1)
  x1 - vector()
  x2 - vector()
  for(i in m){
  x1[i] - i
  x2[i] - i^2
  }
  dat - data.frame(x1,x2)
  But, I have  false result
  dat
x1 x2
  1 1  1

  could some tell me how it is possible to do this?

  Thank you!




-

[[alternative HTML version deleted]]

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

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


Re: [R] regex sulotion for seperating number and string

2008-03-05 Thread sun
Thanks all for the prompt answers!!!  All works perfectly!

up and running! Thanks!

jim holtman [EMAIL PROTECTED] wrote in message 
news:[EMAIL PROTECTED]
 This should do it for you:

 x - c(2564gc, 2367,GH, 2134 JHG)
 x.sep - gsub(([[:digit:]]+)[ ,]*([[:alpha:]]+), \\1 \\2, x)
 # now create separate values
 strsplit(x.sep,  )
 [[1]]
 [1] 2564 gc

 [[2]]
 [1] 2367 GH

 [[3]]
 [1] 2134 JHG




 On 3/5/08, sun [EMAIL PROTECTED] wrote:
 I have strings contain postcode and letters, some seperated with blank, 
 some
 with comma, and some hasn't seperated. eg, 2324gz  2567 HK 3741,BF

 I want to seperate the number and letters into two new variables.

 I know this should be quite basic question, but searched on regex syntax 
 and
 that seems a bit scarey to me, any one can shot me a quick solution on 
 this
 particular question?

 thanks,
 Sun

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



 -- 
 Jim Holtman
 Cincinnati, OH
 +1 513 646 9390

 What is the problem you are trying to solve?

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


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


Re: [R] loop

2008-03-05 Thread Erik Iverson
Why not simply?

m - seq(-1, 1, by = 0.1)
dat - data.frame(m, m^2)

- Erik Iverson

Neuer Arkadasch wrote:
 Hello all,

   I am trying to use 

   m - seq(-1,1,0.1)
   x1 - vector()
   x2 - vector()
   for(i in m){
   x1[i] - i
   x2[i] - i^2
   }
   dat - data.frame(x1,x2)
   But, I have  false result
   dat
 x1 x2
   1 1  1

   could some tell me how it is possible to do this?

   Thank you!


 

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

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


Re: [R] R_alloc with structures with flexible array members

2008-03-05 Thread Jeffrey Horner
Ramon Diaz-Uriarte wrote on 03/05/2008 04:25 AM:
 Dear All,
 
 In a package, I want to use some C code where I am using a structure
 (as the basic element of a linked list) with flexible array members.
 Basically, this is a structure where the last component is an
 incomplete array type  (e.g., Harbison  Steel, C, a reference
 manual, 5th ed., p. 159) such as:
 
 struct Sequence {
   struct Sequence *next;
   int len;
   unsigned int state_count[];
 };
 
 
 To create one such sequence, I allocate storage (following Harbison
 and Steel) in a C program as follows:
 
 struct Sequence *A;
 int n = 4;
 A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int));
 
 
 If I understand correctly, however, it would be better to use R_alloc
 instead of malloc (memory automagically freed on exit and error;
 error-checking). But I do not know how to make the call to R_alloc
 here, since R_alloc allocates n units of size bytes each.
 
 
 I've tried, without success, the following two:
 
 int to_add_for_R_alloc =
 (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int));
 
   A = (struct sequence *) R_alloc(to_add_for_R_alloc + n,
 sizeof(unsigned int));
 
 or even a brute force attempt as:
 
  A = (struct sequence *) R_alloc( 100,  sizeof(struct sequence));
 
 
 but both result in segmentation faults.
 
 
 Should I just keep using malloc (and free at end)?

Hi Ramon,

You should be able to use R_alloc without seg faults, so there's 
something wrong with your code somewhere. R_alloc multiplies its 
arguments together to come up with the total number of bytes to allocate 
then it allocates a raw vector and returns the data portion.

So you can just treat R_alloc similarly to malloc by calling 
R_alloc(1,sizeof(struct Sequence) + n * sizeof(unsigned int)).

Best,

Jeff
-- 
http://biostat.mc.vanderbilt.edu/JeffreyHorner

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


Re: [R] rrp.impute: for what sizes does it work?

2008-03-05 Thread Werner Wernersen
rrp is working!

Sorry, it was my mistake... fiddling around to find
out what the problem is I forgot to re-include the
variables which are to be imputed. It seems like this
case is not caught but the algorithm finishes with the
mentioned error.

Anyway, I am still a little fuzzy about imputation and
browsing the web I found no clear recommendations
regarding this topic. Which R function would you
recommend for imputing categorical variables /
factors?

Many thanks,
  Werner


 Hi,
 
 I have a survey dataset of about 2 observations
 where for 2 factor variables I have about 200
 missing
 values each. I want to impute these using 10
 possibly
 explanatory variables which are a mixture of
 integers
 and factors.
 
 Since I was quite intrigued by the concept of rrp I
 wanted to use it but it takes ages and terminates
 with
 an error. First time it stopped complaining about
 too
 little memory which I increased then from 1.5 to
 2GB. 
 Now it terminates after a long time with this:
 Error in nn[[i]] : subscript out of bounds
 
 Has anybody encountered this problem before, is it
 due
 to my data? Could you recommend another package for
 such imputation? 
 I looked at the aregImpute in Hmisc but I don't
 really
 understand what it is doing.
 
 Thanks a lot,
  Werner



  Lesen Sie Ihre E-Mails jetzt einfach von unterwegs.

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


Re: [R] ipf function in R

2008-03-05 Thread Charles C. Berry
On Wed, 5 Mar 2008, Chandra Shah wrote:

 Hi
 I have a 3 x 2 contingency table:
 10 20
 30 40
 50 60
 I want to update the frequencies to new marginal totals:
 100 130
 40 80 110
 I want to use  the ipf (iterative proportional fitting) function which
 is apparently in the cat package.
 Can somebody please advice me how to input this data and invoke ipf in R
 to obtain an updated contingency table?

I'd use loglin()

newtab -
loglin( rowmarg%o%colmarg/sum(colmarg),
margin=list(1,2),
start=tab, fit=TRUE )$fit


with rowmarg and colmarg set to your updated marginals.

As for inputting the data, if this is all you have, type it in at the 
command line. See

?matrix
?c

and note this:

 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


HTH,

Chuck


 Thanks.
 By the way I am quite new to R.

 -- 


 Dr Chandra Shah
 Senior Research Fellow
 Monash University-ACER Centre for the Economics of Education and Training
 Faculty of Education, Building 6,
 Monash University
 Victoria
 Australia 3800
 Tel. +61 3 9905 2787
 Fax +61 3 9905 9184

 www.education.monash.edu.au/centres/ceet

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


Charles C. Berry(858) 534-2098
 Dept of Family/Preventive Medicine
E mailto:[EMAIL PROTECTED]  UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

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


Re: [R] package for repeated measures ANOVA with various link functions REDUX

2008-03-05 Thread Douglas Bates
On Tue, Mar 4, 2008 at 9:48 PM, John Sorkin [EMAIL PROTECTED] wrote:
 Prof. Bates was correct to point out the lack of specifics in my original 
 posting. I am looking for a package that will allow we to choose among link 
 functions and account for repeated measures in a repeated measures ANOVA.

  My question is what package should I use to facilitate estimating rates of 
 illegal drug use at three centers, and the effect two interventions have on 
 usage. At each center data describing the rate of drug use was obtained once 
 a month. For the first six-months, there was no intervention at any of the 
 three centers. For months seven through 13 intervention one was applied at 
 each of the three centers. For months 14 through 24 intervention two was 
 applied at each center. The question I am trying to answer is did 
 intervention one or two change drug usage at any of the three centers. I am 
 treating center as a repeated measure, i.e. the rate of drug use at month one 
 will be correlated with the rate of drug use at center one at months two, 
 three, etc.

  I have accounted for repeated measures several ways in the past.

  (1) I have used SAS proc MIXED with a REPEATED statement. The REPEATED 
 statement allows for the specification of the within-subject correlation of 
 repeated measures by specifying the structure of the within-subject 
 variance-covariance matrix of the repeated measures. The matrix is block 
 diagonal with one block for each subject.

But does such a structure extend to models with binary or count
responses?  You have mentioned that you want to use an arbitrary link
function such as quasibinomial.  What I understand the effect of the
REPEATED statement to be is to specify a parameterized form of the
marginal variance-covariance matrix of the responses.  If the response
variable has a multivariate normal distribution it is possible to
independently specify the mean (determined by the fixed-effects
parameters) and the marginal variance-covariance.

However, in the case of generalized linear models the mean response is
determined by a linear predictor and a link function while the
variance-covariance of the response is determined by prior weights and
a variance function.  The same is true for generalized linear mixed
models except that this description applies to the conditional
distribution of the response given the random effects.  The link and
the variance functions must agree so, for example, using a logit or
probit link which restricts the value of mu to the interval [0,1]
would imply a variance function (up to prior weights) of mu(1-mu).  At
least I think so - others may feel that it is possible to specify an
arbitrary variance function but I don't see how that can make sense.
To me the whole point of generalized linear models is to transform the
linear predictor to the desired range for the mean and to take into
account what this implies about the variance.

Even if you feel that it is possible to relax the ties between the
link function and the variance function I don't see how it would be
possible to specify an arbitrary structure for the marginal
variance-covariance of the response.  If you say that the marginal
variance-covariance must have a block-wise compound symmetry structure
but you are going to restrict the mean to the range [0,1] I think you
have painted yourself into a corner.  I don't think it is possible to
specify a mean on a restricted range and separately specify an
arbitrary variance-covariance structure.  In particular, when the mean
is on the range [0,1] then you better have the variance going to zero
as the mean goes to 0 or to 1.  You can't arbitrarily say that the
variance within a block must be constant, regardless of the values of
the means in those blocks.


  (2) I have used SAS proc GENMOD which uses GEE to adjust the parameter 
 estimates and their standard errors for the fact that a repeated measurements 
 of a parameter are obtained from a given subjects.

  Is there any package in R that will allow me to perform a repeated measures 
 ANOVA with a selection of link functions that will allow me to account for 
 repeated measures by either specifying the correlation structure of the 
 repeated measures from a subject a la SAS proc mixed or by adjusting the 
 parameter estimates using GEE a la proc GENMOD? Perhaps R has a package that 
 accounts for repeated measures in some other manner.

  Thank you,
  John Sorkin



  John Sorkin M.D., Ph.D.
  Chief, Biostatistics and Informatics
  University of Maryland School of Medicine Division of Gerontology
  Baltimore VA Medical Center
  10 North Greene Street
  GRECC (BT/18/GR)
  Baltimore, MD 21201-1524
  (Phone) 410-605-7119
  (Fax) 410-605-7913 (Please call phone number above prior to faxing)

   Douglas Bates [EMAIL PROTECTED] 3/4/2008 5:13 PM 
  On Tue, Mar 4, 2008 at 10:52 AM, John Sorkin
  [EMAIL PROTECTED] wrote:
   R 2.6.0
Windows XP

At the risk of raising the ire of the R gods . . 

[R] Correlation matrix one side with significance

2008-03-05 Thread Martin Kaffanke
Hi there!

In my case, 

cor(d[1:20])

makes me a good correlation matrix.

Now I'd like to have it one sided, means only the left bottom side to be
printed (the others are the same) and I'd like to have * where the
p-value is lower than 0.05 and ** lower than 0.01.

How can I do this?

And another thing: Is there a way to output that table as a latex table?

Thanks,
Martin

-- 
Ihr Partner für Webdesign, Webapplikationen und Webspace.
http://www.roomandspace.com/
Martin Kaffanke +43 650 4514224


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


Re: [R] Asking, are simple effects different from 0

2008-03-05 Thread Chuck Cleland
On 3/5/2008 10:09 AM, jebyrnes wrote:
 Huh.  Very interesting.  I haven't really worked with manipulating contrast
 matrices before, save to do a prior contrasts.  Could you explain the matrix
 you laid out just a bit more so that I can generalize it to my case?  

   Each column corresponds to one of the coefficients in the model, and 
each row specifies a particular contrast.  The numbers in the matrix 
indicate how the model coefficients are combined to indicate a 
particular difference in means.
   For example, the first row indicates that the third coefficient 
(woolB) is multiplied by -1.  The baseline categories are A and L for 
the wool and tension factors, so the woolB effect in fm is the simple 
effect of B vs. A in the baseline category of the tension factor. 
Multiplying this coefficient by -1 produces an A vs. B comparison in the 
baseline category of the tension factor.
   When I want to check that contrasts are as intended, I use contrast() 
in the contrast package (by Steve Weston, Jed Wing,  Max Kuhn).  That 
function allows you to specify factor levels by name to construct the 
contrast.  For example:

library(contrast)

# M vs. H at B

contrast(fm, a=list(tension = M, wool = B),
  b=list(tension = H, wool = B))

lm model parameter contrast

  Contrast S.E.  LowerUppert df Pr(|t|)
10 5.157299 -0.3694453 20.36945 1.94 48   0.0584

   It also allows you to print the design matrix for a contrast:

contrast(fm, a=list(tension = M, wool = B),
  b=list(tension = H, wool = B))$X

   (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB
1   01   -1 0  1 -1

 Chuck Cleland wrote:

One approach would be to use glht() in the multcomp package.  You 
 need to work out how to formulate the matrix of coefficients that give 
 the desired contrasts.  Here is an example using the warpbreaks data
 frame:

 fm - lm(breaks ~ tension*wool, data=warpbreaks)

 # names(coef(fm))
 # (Intercept) tensionM tensionH woolB tensionM:woolB tensionH:woolB

 cm - rbind(
 A vs. B at L = c(0, 0, 0,-1, 0, 0),
 A vs. B at M = c(0, 0, 0,-1,-1, 0),
 A vs. B at H = c(0, 0, 0,-1, 0,-1),
 M vs. L at A = c(0, 1, 0, 0, 0, 0),
 M vs. H at A = c(0, 1,-1, 0, 0, 0),
 L vs. H at A = c(0, 0,-1, 0, 0, 0),
 M vs. L at B = c(0, 1, 0, 0, 1, 0),
 M vs. H at B = c(0, 1,-1, 0, 1,-1),
 L vs. H at B = c(0, 0,-1, 0, 0,-1))

 library(multcomp)

 summary(glht(fm, linfct = cm), test = adjusted(type=none))

   Simultaneous Tests for General Linear Hypotheses

 Fit: lm(formula = breaks ~ tension * wool, data = warpbreaks)

 Linear Hypotheses:
Estimate Std. Error t value  p value
 A vs. B at L == 0  16. 5.1573   3.167 0.002677 **
 A vs. B at M == 0  -4.7778 5.1573  -0.926 0.358867
 A vs. B at H == 0   5.7778 5.1573   1.120 0.268156
 M vs. L at A == 0 -20.5556 5.1573  -3.986 0.000228 ***
 M vs. H at A == 0  -0.5556 5.1573  -0.108 0.914665
 L vs. H at A == 0  20. 5.1573   3.878 0.000320 ***
 M vs. L at B == 0   0.5556 5.1573   0.108 0.914665
 M vs. H at B == 0  10. 5.1573   1.939 0.058392 .
 L vs. H at B == 0   9. 5.1573   1.831 0.073270 .
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 (Adjusted p values reported -- none method) 

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


Re: [R] inheritence in S4

2008-03-05 Thread cgenolin
Well well well...

To summarize : let assume that A is a class (slot x) and C is a class 
containing A (A and slot y) - as(c,A) calls new(A). So new(A) HAS 
TO works, you can not decide to forbid empty object (unless you define 
setAs(C,A) ?)
- In addition, any test that you would like to set in initialize or 
validity have to first check is some field are empty (because 
'if([EMAIL PROTECTED] 0)' will fail if x=numerical(0))
- When you call new(C), the neither new(A) nor intialize(A) are 
called. ()

So, all the security test you write in initialize A, you have to 
rewrite them on C ?

I start to undestand why not that much people use S4...

Christophe

 [EMAIL PROTECTED] writes:

 Hi Martin, thanks for your answer

 But a couple
 of other quick points. I would have written

 setMethod(initialize, A,
  function(.Object, ..., xValue=numeric(0)){
  callNextMethod(.Object, ..., x=xValue)
  })


 I am not that much familiar with S3... In our way of writing this
 method, 'initialize' for 'A' will call the next method for A, which is
 initialise' for 'numeric', is that right ?

 from below,

 setClass(A, representation(x=numeric))

 The 'next' method refers to arguments in the signature of the generic,
 i.e., the 'next' method relevant to .Object. For A, the next method
 is the default initialize method, which will use all named arguments
 to fill the slots of .Object (and perhaps call validObject, among
 other things). You can see the code with

 showMethods(initialize, classes=ANY, includeDef=TRUE)

 S3 is not relevant here.

 And finally, the position of 'xValue' and 'yValue' means that
 the arugment has to be named, e.g., new(B, yValue=12). This seems a
 little awkward at first, but seems like a best practice

 I agree with you. But I do not like the use of ... , it lets us to
 make many mistake like in :

 print(3.5165,digitts=2)

 There is a typo in digitts, but since print works with ... , R does
 not detect the mistake.
 So I agree with the importance of naming argument, I always do it but
 without 'forcing' it

 If this were initialize, and you had provided an invalid named
 argument, the default method would have failed because there is no
 slot of that name.

 And finally, in Jim's thread I mention using a constructor. So in
 practice for a case like the above I would not define any initialize
 methods,

 Interesting. Will you still define a validity method or not even ?

 For a simple case like this there is no extra validity to check beyond
 that ensured by S4 (e.g., slots of the correct type). If there were
 constraints, e.g., only positive values, or length one vectors, then
 I would define a validity function.

 Martin


 Christophe



 Christophe Genolini [EMAIL PROTECTED] writes:

 Thanks Martin

 Well it works except that as seems to not like the initialize
 method : the following code (that is the same than yours with some
 initialize for A B and C) does not compile. It seems that as(c,A)
 does not work if we definie a initialize for A...

 --- 8 --
 setClass(A, representation(x=numeric))
 setMethod(initialize,A,function(.Object,value)[EMAIL PROTECTED] -
 value;return(.Object)})
 a - new(A,4)

 setClass(B, representation(y=numeric))
 setMethod(initialize,B,function(.Object,value)[EMAIL PROTECTED] -
 value;return(.Object)})
 b - new(B,5)

 setClass(C, contains=c(A, B))
 setMethod(initialize,C,function(.Object,valueA, valueB){
 [EMAIL PROTECTED] - valueA
 [EMAIL PROTECTED] - valueB
 return(.Object)
 })
 c - new(C,valueA=10,valueB=12)

 setMethod(show, A, function(object) cat(A\n))
 setMethod(show, B, function(object) cat(B\n))
 setMethod(show, C, function(object) {
 callGeneric(as(object, A))
 callGeneric(as(object, B))
 cat(C\n)
 })
 c
 --- 8 

 Is there something wrong with the use of 'as' between class and father
 class?

 Christophe
 Hi Christophe --

 I don't know whether there's a particularly elegant way. This works

 setClass(A, representation(x=numeric))
 setClass(B, representation(y=numeric))
 setClass(C, contains=c(A, B))

 setMethod(show, A, function(object) cat(A\n))
 setMethod(show, B, function(object) cat(B\n))
 setMethod(show, C, function(object) {
 callGeneric(as(object, A))
 callGeneric(as(object, B))
 cat(C\n)
 })


 new(C)

 A
 B
 C

 but obviously involves the developer in making explicit decisions
 about method dispatch when there is multiple inheritance.

 Martin

 [EMAIL PROTECTED] writes:


 Hi the list

 I define a class A (slot a and b), a class C (slot c and d) and a
 class E that inherit from A and B.
 I define print(A) and print(B). For print(C), I would like to use
 both of them, but I do not see how...

 Thanks for your help...

 Christophe

 
 Ce message a ete envoye par IMP, grace a l'Universite Paris 10 Nanterre

 __
 R-help@r-project.org mailing list
 

Re: [R] Reversed but positive axis in trellis plots?

2008-03-05 Thread Deepayan Sarkar
On 3/5/08, Fredrik Karlsson [EMAIL PROTECTED] wrote:
 Hi,

  In my discpipline, it is common to plot one acoustic property on a
  positive scale but from top to bottom on the ordinate and the same for
  another measurement on the abscissa.
  So, the origin of the plot is on the top right of the plot, with
  increasing values to the left /down. This is to highlight the
  correlation between the acoustic measurement and the position of the
  forming structure, for instance when teaching it to students.

  The grouping ability of the trellis plot is quite handy whan plotting
  many instances of the same thing, so I was wondering if it is possible
  to make a trellis xyplot behave this way?
  Converting all values to negative and changing labels to the negative
  of the negative seems one solution to the reverseness of the axes, but
  how do I change the position? Is it possible?

You can always specify explicit limits that are reversed, e.g.,
ylim=c(100, 0). If you want this as a general behaviour but don't know
your data's range beforehand, one option might have been

prepanel = function(x, y, ...) {
list(ylim = rev(range(y)))
}

Unfortunately, this doesn't work when relation=same, as the step of
combining the per-panel limits to obtain a common range disregards the
order.

-Deepayan

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


Re: [R] loop

2008-03-05 Thread Neuer Arkadasch
Thank you Yinghai, that's what I need :-)!

Yinghai Deng [EMAIL PROTECTED] schrieb: m - seq(-1,1,0.1)
  x1 - c()
  x2 - c()
  for(i in 1:length(m)){
  x1[i] - m[i]
  x2[i] - m[i]^2
  }
  dat - data.frame(x1,x2)

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Behalf Of Neuer Arkadasch
Sent: March 5, 2008 10:00 AM
To: [EMAIL PROTECTED]
Subject: [R] loop


Hello all,

  I am trying to use

  m - seq(-1,1,0.1)
  x1 - vector()
  x2 - vector()
  for(i in m){
  x1[i] - i
  x2[i] - i^2
  }
  dat - data.frame(x1,x2)
  But, I have  false result
  dat
x1 x2
  1 1  1

  could some tell me how it is possible to do this?

  Thank you!




-

 [[alternative HTML version deleted]]

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



   
-
  Lesen Sie Ihre E-Mails jetzt einfach von unterwegs.. 
[[alternative HTML version deleted]]

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


[R] density plot

2008-03-05 Thread Tom Hart
Hi,

I'm trying to create a density plot which I used to do in geneplotter
using the following code. Unfortunately I can't find the combination of
R release and geneplotter that works. 

Can anyone suggest a fix or an alternative to smoothScatter that will
plot depth of one dive vs depth of the next as density?

Many thanks,

Tom

 

library(geneplotter)

require(RColorBrewer)

layout(matrix(1:4, ncol=2, byrow=TRUE))

op - par(mar=rep(2,4))

 

  x   - cbind(depthcurrent1,depthcurrent2)

  smoothScatter(x,nrpoints=0, colramp=colorRampPalette(c(white,
black)), ylab=Depth of Next Dive, xlab=Depth of Current Dive)

 

  colors  - densCols(x)

  plot(x, col=black, pch=20, ylab=Depth of Next Dive, xlab=Depth of
Current Dive)

 



The Zoological Society of London is incorporated by Royal Charter
Principal Office England. Company Number RC000749
Registered address: 
Regent's Park, London, England NW1 4RY
Registered Charity in England and Wales no. 208728 

_
This e-mail has been sent in confidence to the named add...{{dropped:20}}

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


Re: [R] CROSSOVER TRIALS IN R (Binary Outcomes)

2008-03-05 Thread Charles C. Berry
On Wed, 5 Mar 2008, Boikanyo Makubate wrote:

 I will like to analyse a binary cross over design using the random
 effects model. The probability of success is assumed to be logistic.
 Suppose as an example, we have 4 subjects undergoing a crossover design,
 where the outcome is either success or failure. The first two subjects
 receive treatment A first followed by treatment B. The remaining two
 subjects receive treatments in the reverse order. The outcomes for the
 subjects is the sequence AB are as follows: (0,1) and (0,0). While the
 outcomes for the subjects in the sequence BA are (1,1) and (1,0).

 How can i analyse this using R. I have done the problem with PROC
 NLMIXED in SAS, I simply want to compare the results from SAS with those
 from R. Please help. R-Codes will be highly appreciated.


You need to be clear about the statistical model.

In crossover trials, there is usually a subject effect that one conditions 
out of the likelihood (often implicitly).

In this case, clogit() in the survival package would be suitable for such 
an approach. Something like

res - clogit( outcomes ~ rx = strata( subject ) )

However, the example you give is degenerate in the same way that an 
ordinary logistic regression is when the responses are linearly separable 
in the regressor space.

So, a well crafted program will tell you that it cannot find an answer 
(and clogit does this) for your example data.

If you knew all of this and had some other statistical model in mind 
(such as one that models the subject effects rather than conditioning 
on them), you should repost telling us what model you wanted to fit.

HTH,

Chuck


 Boikanyo.

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


Charles C. Berry(858) 534-2098
 Dept of Family/Preventive Medicine
E mailto:[EMAIL PROTECTED]  UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

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


Re: [R] Correlation matrix one side with significance

2008-03-05 Thread Henrique Dallazuanna
Try this:

On 05/03/2008, Martin Kaffanke [EMAIL PROTECTED] wrote:
 Hi there!

  In my case,

  cor(d[1:20])

  makes me a good correlation matrix.

  Now I'd like to have it one sided, means only the left bottom side to be
  printed (the others are the same) and I'd like to have * where the
  p-value is lower than 0.05 and ** lower than 0.01.

  How can I do this?

d - matrix(rexp(16, 2), 4)
corr - cor(d)
sign - symnum(cor(d), cutpoints=c(0.05, 0.01), corr = T,
symbols=c(***, **, *), abbr=T, diag=F)

noquote(mapply(function(x, y)paste(x, format(y, dig=3), sep=''),
as.data.frame(unclass(sign)), as.data.frame(corr)))



  And another thing: Is there a way to output that table as a latex table?

See ?latex function in Hmisc package and also xtable package

  Thanks,
  Martin


  --
  Ihr Partner für Webdesign, Webapplikationen und Webspace.
  http://www.roomandspace.com/
  Martin Kaffanke +43 650 4514224

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





-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40 S 49° 16' 22 O

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


Re: [R] vector manipulations

2008-03-05 Thread Pete Dorothy
Thank you everybody.

Phil, your expand.grid works very nicely and I will use it for
non-vectorized functions.

Yet I am a bit confused about vectorization. For me it is synonymous of
no loop. :-(

I wrote a toy example (with a function which is not my log-likelihood).

FIRST PART

nir=1:10
logl=function(x,y,nir) sum(log(x*nir+y))

x=seq(0.1,0.3,by=10^(-1))
y=seq(0.1,0.3,by=10^(-1))
z=outer(x,y,logl,nir=nir)

This does not work. Can you explain me why it is not vectorised ?

SECOND PART

nir=1:10
logl2=function(x,y,nir) {
a=0
for (i in 1:10) {
a=a+log(x*nir[i]+y)
}
return(a)
}

x=seq(0.1,0.3,by=10^(-1))
y=seq(0.1,0.3,by=10^(-1))
z2=outer(x,y,logl2,nir=nir)

This seems to work, though the function does not seem to be vectorized.

I am sorry for being such a noob. I'm ok in maths but I am bad at
programming. I bought a book on R (Introductory Statistics with R by
Dalgaard) ** on Amazon last week . I will read it when I receive it. Do you
know other good books ?

2008/3/5, [EMAIL PROTECTED] [EMAIL PROTECTED]:

 No problems with it working.  The main problem I have observed is
 unrealistic expectations.  People write an *essentially* non-vectorized
 function and expect Vectorize to produce a version of it which will
 out-perform explicit loops every time.  No magic bullets in this game.

 Bill.



 Bill Venables
 CSIRO Laboratories
 PO Box 120, Cleveland, 4163
 AUSTRALIA
 Office Phone (email preferred): +61 7 3826 7251
 Fax (if absolutely necessary):  +61 7 3826 7304
 Mobile: +61 4 8819 4402
 Home Phone: +61 7 3286 7700
 mailto:[EMAIL PROTECTED]
 http://www.cmis.csiro.au/bill.venables/

 -Original Message-

 From: Duncan Murdoch [mailto:[EMAIL PROTECTED]
 Sent: Wednesday, 5 March 2008 9:36 AM
 To: Venables, Bill (CMIS, Cleveland)
 Cc: r-help@r-project.org
 Subject: Re: [R] vector manipulations

 On 3/4/2008 5:41 PM, [EMAIL PROTECTED] wrote:
  Your problem is that your function log1( , ) is not vectorized with
  respect to its arguments.  For a function to work in outer(...) it
 must
  accept vectors for its first two arguments and it must produce a
  parallel vector of responses.
 
  To quote the help information for outer:
 
  FUN is called with these two extended vectors as arguments.
 Therefore,
  it must be a vectorized function (or the name of one), expecting at
  least two arguments.
 
  Sometimes Vectorize can be used to make a non-vectorized function into
 a
  vectorized one, but the results are not always entirely satisfactory
 in
  my experience.

 What problems have you seen?

 Duncan Murdoch

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


[[alternative HTML version deleted]]

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


Re: [R] R_alloc with structures with flexible array members

2008-03-05 Thread Ramon Diaz-Uriarte
Dear Jeff,

Thanks for the suggestion. However, something is still not working.
This is a simple example:

***  start C 
#include R.h

struct Sequence {
  int len;
  unsigned int state_count[];
};


int main(void) {

  struct Sequence *A;
  int n = 4;

  // First line segfaults. Second doesn't
A = (struct Sequence *) R_alloc(1,  sizeof(struct Sequence) + n *
sizeof(unsigned int));
  //  A = malloc(sizeof(struct Sequence) + n * sizeof(unsigned int));

  return(0);
}

***  end C  **


I then do
gcc -std=gnu99 -Wall -I/usr/share/R/include -I/usr/share/R/include
-L/usr/lib/R/lib -lR ex7.c

and the ./a.out segfaults when I use R_alloc (not with malloc).


Best,

R.

On Wed, Mar 5, 2008 at 5:23 PM, Jeffrey Horner
[EMAIL PROTECTED] wrote:
 Ramon Diaz-Uriarte wrote on 03/05/2008 04:25 AM:


  Dear All,
  
   In a package, I want to use some C code where I am using a structure
   (as the basic element of a linked list) with flexible array members.
   Basically, this is a structure where the last component is an
   incomplete array type  (e.g., Harbison  Steel, C, a reference
   manual, 5th ed., p. 159) such as:
  
   struct Sequence {
 struct Sequence *next;
 int len;
 unsigned int state_count[];
   };
  
  
   To create one such sequence, I allocate storage (following Harbison
   and Steel) in a C program as follows:
  
   struct Sequence *A;
   int n = 4;
   A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int));
  
  
   If I understand correctly, however, it would be better to use R_alloc
   instead of malloc (memory automagically freed on exit and error;
   error-checking). But I do not know how to make the call to R_alloc
   here, since R_alloc allocates n units of size bytes each.
  
  
   I've tried, without success, the following two:
  
   int to_add_for_R_alloc =
   (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int));
  
 A = (struct sequence *) R_alloc(to_add_for_R_alloc + n,
   sizeof(unsigned int));
  
   or even a brute force attempt as:
  
A = (struct sequence *) R_alloc( 100,  sizeof(struct sequence));
  
  
   but both result in segmentation faults.
  
  
   Should I just keep using malloc (and free at end)?

  Hi Ramon,

  You should be able to use R_alloc without seg faults, so there's
  something wrong with your code somewhere. R_alloc multiplies its
  arguments together to come up with the total number of bytes to allocate
  then it allocates a raw vector and returns the data portion.

  So you can just treat R_alloc similarly to malloc by calling
  R_alloc(1,sizeof(struct Sequence) + n * sizeof(unsigned int)).

  Best,

  Jeff
  --
  http://biostat.mc.vanderbilt.edu/JeffreyHorner




-- 
Ramon Diaz-Uriarte
Statistical Computing Team
Structural Biology and Biocomputing Programme
Spanish National Cancer Centre (CNIO)
http://ligarto.org/rdiaz

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


[R] Replace values in data.frame conditional on another data.frame

2008-03-05 Thread Carson Farmer
Dear List,

I am looking for an efficient method for replacing values in a
data.frame conditional on the values of a separate data.frame. Here is
my scenario:
I have a data.frame (A) with say 1000 columns, and 365 rows. Each cell
in the data.frame has either valid value, or NA. I have an additional
data.frame (B) with the same number of rows and columns, with valid
values in all cells. What I would like to do, is replace the cells in B
with NA if and only if the corresponding cell in A is NA.
I have search extensively for a method to do this, and I'm sure that in
the end the solution will be embarrassingly simple, however, any help is
greatly appreciated!

Cheers,

Carson

Apologies if this has been posted twice, we are currently experiencing 
server problems...

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


Re: [R] Replace values in data.frame conditional on another data.frame

2008-03-05 Thread John Kane
Try 
bb[is.na(aa)] - NA

It may be simple but it is not necessarily obvious :)
--- Carson Farmer [EMAIL PROTECTED] wrote:

 Dear List,
 
 I am looking for an efficient method for replacing
 values in a
 data.frame conditional on the values of a separate
 data.frame. Here is
 my scenario:
 I have a data.frame (A) with say 1000 columns, and
 365 rows. Each cell
 in the data.frame has either valid value, or NA. I
 have an additional
 data.frame (B) with the same number of rows and
 columns, with valid
 values in all cells. What I would like to do, is
 replace the cells in B
 with NA if and only if the corresponding cell in A
 is NA.
 I have search extensively for a method to do this,
 and I'm sure that in
 the end the solution will be embarrassingly simple,
 however, any help is
 greatly appreciated!
 
 Cheers,
 
 Carson
 
 Apologies if this has been posted twice, we are
 currently experiencing 
 server problems...
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained,
 reproducible code.


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


[R] New data source - now how do we build an R interface?

2008-03-05 Thread Ajay Shah
Folks,

A nice new data resource has come up -- http://data.un.org/
I thought it would be wonderful to setup an R function like
tseries::get.hist.quote() which would be able to pull in some or all
of this data.

I walked around a bit of it and I'm not able to map the resources to
predictable URLs which can then be wget. There's some javascript going
on that I'm not understanding. Perhaps someone on the mailing list can
think about this?

-- 
Ajay Shah  http://www.mayin.org/ajayshah  
[EMAIL PROTECTED] http://ajayshahblog.blogspot.com
*(:-? - wizard who doesn't know the answer.

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


Re: [R] loop

2008-03-05 Thread John Kane
m - seq(-1,1,0.1)
  x1 - vector()
  x2 - vector()
# the loop statement was incorrect. 
  for(i in 1:length(m)){   
  x1[i] - m[i]
  x2[i] - m[i]^2
  }
  dat - data.frame(x1,x2)

# But why not something like this? There is no need
for a loop.

x1 - seq(-1,1,0.1)
mdat - data.frame(x1, x2=x1^2)
--- Neuer Arkadasch [EMAIL PROTECTED] wrote:

 Hello all,

   I am trying to use 

   m - seq(-1,1,0.1)
   x1 - vector()
   x2 - vector()
   for(i in m){
   x1[i] - i
   x2[i] - i^2
   }
   dat - data.frame(x1,x2)
   But, I have  false result
   dat
 x1 x2
   1 1  1

   could some tell me how it is possible to do this?

   Thank you!


 

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


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


Re: [R] Asking, are simple effects different from 0

2008-03-05 Thread jebyrnes

Ah.  I see.  So, if I want to test to see whether each simple effect is
different from 0, I would do something like the following:


cm2 - rbind( 
A:L =  c(1, 0, 0, 0, 0, 0), 
A:M = c(1, 1, 0, 0, 0, 0), 
A:H = c(1, 0, 1, 0, 0, 0), 
B:L =   c(1, 0, 0, 1, 0, 0), 
B:M = c(1, 1, 0, 1, 1, 0), 
B:H = c(1, 0, 1, 1, 0, 1))

summary(glht(fm, linfct = cm2), test = adjusted(type=none)) 

Correct? What is the df on those t-tests then?  Is it 48?

Interestingly, I find this produces results no different than

fm2-lm(breaks ~ tension:wool+0, data=warpbreaks) 
summary(fm2)

Also, here, it would seem each t-test was done with the full 48df.  Hrm.


Chuck Cleland wrote:
 
 
Each column corresponds to one of the coefficients in the model, and 
 each row specifies a particular contrast.  The numbers in the matrix 
 indicate how the model coefficients are combined to indicate a 
 particular difference in means.
For example, the first row indicates that the third coefficient 
 (woolB) is multiplied by -1.  The baseline categories are A and L for 
 the wool and tension factors, so the woolB effect in fm is the simple 
 effect of B vs. A in the baseline category of the tension factor. 
 Multiplying this coefficient by -1 produces an A vs. B comparison in the 
 baseline category of the tension factor.
  
 
 

-- 
View this message in context: 
http://www.nabble.com/Asking%2C-are-simple-effects-different-from-0-tp15835552p15855630.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] box-constrained

2008-03-05 Thread Gustave Lefou
Hello everybody,

I have a question about box-constrained optimization. I've done some
research and I found that optim could do that. Are there other ways in R ?

Is the following correct if I have a function f of two parameters belonging
for example to [0,1] and [0,Infinity] ?
optim(par=param, fn=f, method=L-BFGS-B, lower=c(0,0), upper=c(1,Inf))

My other question is whether it is possible to add the derivatives of my
function (like in nlm) and whether it is better to add them ?

If there is no need to add the derivatives, then I guess I could wish to
optimize the likelihood directly rather than to optimize the
log-likelihood... Indeed one aspect of the log-likelihood is to make the
derivatives tractable (in iid cases). Do you agree ?

Thank you !
Gustave

[[alternative HTML version deleted]]

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


Re: [R] legend for several graphics

2008-03-05 Thread John Kane
?mtitle should do it.
--- Georg Otto [EMAIL PROTECTED] wrote:

 Hi,
 
 I am trying to generate a figure of 9 plots that are
 contained in one
 device by using
 
 par(mfrow = c(3,3,))
 
 I would like to have 1 common legend for all 9 plots
 somewhere outside
 of the plotting area (as opposed to one legend
 inside each of the 9
 plots, which the function legend() seems to generate
 by default).
 
 Any hint how to do this?
 
 Best,
 
 Georg
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained,
 reproducible code.


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


Re: [R] simulation study using R

2008-03-05 Thread Davood Tofighi
Thanks to All,

The comments were very helpful; however, the the simulation is running very
slow. I reduced the number of loops (conditions) so I have 36 loops, and the
data-generation occurs 1000 times within each loop. At the end of each 1000
reps, I saved the summary (e.g., mean) of the reps to a single row vector,
and then saved it in a fileshare database. When the simulation is
finished,  I rbind() the 36 rows of the database objects into a final
simulation result.

Thanks again,
Davood Tofighi

On Wed, Mar 5, 2008 at 8:44 AM, Paul Gilbert [EMAIL PROTECTED]
wrote:

 Davood Tofighi wrote:
  Thanks for your reply. For each condition, I will have a matrix or data
  frames of 1000 rows and 4 columns. I also have a total of 64 conditions
 for
  now. So, in total, I will have 64 matrices or data frames of 1000 rows
 and 4
  columns. The format of data I would like to store would be data frames
 or
  matrices. I also would like to store the data for later use,
 I generally find it is better to store the seed and other data you need
 to reproduce the experiment, rather than trying to store the data (see,
 for example, package setRNG). If you save only the summary statistics
 you need, then you can usually do it in memory. (Be sure to assign
 variables for the statistics to their full size first and then populate
 them, rather than extending them at each step.) If you write things to
 files then it will slow down your simulation a lot. In fact, in most
 cases it will be quicker to re-run the experiment than it will be to
 read the data from disk.

 Paul
  e.g., a plot of
  the empirical distribution of the chi^2, or to compute the power of
 Chi^2
  across 1000 reps for each condition.
 
  On Mon, Mar 3, 2008 at 7:03 PM, jim holtman [EMAIL PROTECTED] wrote:
 
 
  What is the format of the data you are storing (single value,
  multivalued vector, matrix, dataframe, ...)?  This will help formulate
  a solution.  What do you plan to do with the data?  Are you going to
  do further analysis, write it to flat files, store it in a data base,
  etc.?  How big are the data objects you are manipulating?
 
  On Mon, Mar 3, 2008 at 7:05 PM, Davood Tofighi [EMAIL PROTECTED]
 wrote:
 
  Dear All,
 
  I am running a Monte Carlo simulation study and have some questions on
 
  how
 
  to manage data storage efficiently at the end of each 1000 replication
 
  loop.
 
  I have three conditions coded using the FOR {} loops and a FOR loop
 that
  generates data for each condition, performs analysis, and computes a
  statistic 1000 times. Therefore, for each condition, I will have 1000
  statistic values. My question is what's the best way to store the 1000
  statistic for each condition. Any suggestion on how to manage such
  simulation studies is greatly appreciated.
  Thanks,
 
  --
  Davood Tofighi
  Department of Psychology
  Arizona State University
 
 [[alternative HTML version deleted]]
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
 
  http://www.R-project.org/posting-guide.html
 
  and provide commented, minimal, self-contained, reproducible code.
 
 
 
  --
  Jim Holtman
  Cincinnati, OH
  +1 513 646 9390
 
  What is the problem you are trying to solve?
 
 
 
 
 
 

 

 La version française suit le texte anglais.


 

 This email may contain privileged and/or confidential information, and the
 Bank of
 Canada does not waive any related rights. Any distribution, use, or
 copying of this
 email or the information it contains by other than the intended recipient
 is
 unauthorized. If you received this email in error please delete it
 immediately from
 your system and notify the sender promptly by email that you have done so.


 

 Le présent courriel peut contenir de l'information privilégiée ou
 confidentielle.
 La Banque du Canada ne renonce pas aux droits qui s'y rapportent. Toute
 diffusion,
 utilisation ou copie de ce courriel ou des renseignements qu'il contient
 par une
 personne autre que le ou les destinataires désignés est interdite. Si vous
 recevez
 ce courriel par erreur, veuillez le supprimer immédiatement et envoyer
 sans délai à
 l'expéditeur un message électronique pour l'aviser que vous avez éliminé
 de votre
 ordinateur toute copie du courriel reçu.




-- 
Davood Tofighi
Department of Psychology
Arizona State University
P.O. BOX 871104
Tempe, AZ 85287-1104
Tel.:480-727-7884

[[alternative HTML version deleted]]

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

Re: [R] Asking, are simple effects different from 0

2008-03-05 Thread Chuck Cleland
On 3/5/2008 1:32 PM, jebyrnes wrote:
 Ah.  I see.  So, if I want to test to see whether each simple effect is
 different from 0, I would do something like the following:
 
 cm2 - rbind( 
 A:L =  c(1, 0, 0, 0, 0, 0), 
 A:M = c(1, 1, 0, 0, 0, 0), 
 A:H = c(1, 0, 1, 0, 0, 0), 
 B:L =   c(1, 0, 0, 1, 0, 0), 
 B:M = c(1, 1, 0, 1, 1, 0), 
 B:H = c(1, 0, 1, 1, 0, 1))

   That does not corresponds to what I think of as the simple effects. 
That specifies the six cell means, but it does not *compare* any cell 
means.  I think of a simple effect as the effect of one factor at a 
specific level of some other factor.

 summary(glht(fm, linfct = cm2), test = adjusted(type=none)) 
 
 Correct? What is the df on those t-tests then?  Is it 48?

   Yes, df = 48 for each contrast.

 Interestingly, I find this produces results no different than
 
 fm2-lm(breaks ~ tension:wool+0, data=warpbreaks) 
 summary(fm2)

   Yes, but those are not what I would call the simple effects.  Those 
are essentially one-sample t-tests for each of the 6 cell means.

 Also, here, it would seem each t-test was done with the full 48df.  Hrm.

   The df are based on the whole model, not the 9 observations in one cell.

 Chuck Cleland wrote:

Each column corresponds to one of the coefficients in the model, and 
 each row specifies a particular contrast.  The numbers in the matrix 
 indicate how the model coefficients are combined to indicate a 
 particular difference in means.
For example, the first row indicates that the third coefficient 
 (woolB) is multiplied by -1.  The baseline categories are A and L for 
 the wool and tension factors, so the woolB effect in fm is the simple 
 effect of B vs. A in the baseline category of the tension factor. 
 Multiplying this coefficient by -1 produces an A vs. B comparison in the 
 baseline category of the tension factor. 

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


[R] N-way ANOVA

2008-03-05 Thread Paul Smith
Dear All,

Can R perform n-way ANOVA, i.e., with 3 or more factors?

Thanks in advance,

Paul

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


Re: [R] legend for several graphics

2008-03-05 Thread Gavin Simpson
On Wed, 2008-03-05 at 15:28 +0100, Georg Otto wrote:
 Hi,
 
 I am trying to generate a figure of 9 plots that are contained in one
 device by using
 
 par(mfrow = c(3,3,))
 
 I would like to have 1 common legend for all 9 plots somewhere outside
 of the plotting area (as opposed to one legend inside each of the 9
 plots, which the function legend() seems to generate by default).
 
 Any hint how to do this?

Here's one way:

op - par(mfrow = c(3,3), ## split region
  oma = c(5,0,4,0) + 0.1, ## create outer margin
  mar = c(5,4,2,2) + 0.1) ## shrink some margins
plot(1:10, main = a, pch = 1:2, col= 1:2)
plot(1:10, main = b, pch = 1:2, col= 1:2)
plot(1:10, main = c, pch = 1:2, col= 1:2)
plot(1:10, main = d, pch = 1:2, col= 1:2)
plot(1:10, main = e, pch = 1:2, col= 1:2)
plot(1:10, main = f, pch = 1:2, col= 1:2)
plot(1:10, main = g, pch = 1:2, col= 1:2)
plot(1:10, main = h, pch = 1:2, col= 1:2)
plot(1:10, main = i, pch = 1:2, col= 1:2)
## title
mtext(My Plots, side = 3, outer = TRUE, font = 2, line = 1, cex = 1.2)
## draw legend
legend(-12.5, -6, legend = c(Type 1, Type 2), pch = 1:2, col = 1:2,
ncol = 2)
par(op)

I had to fiddle by hand with the legend x and y locations to get it
roughly centred. There has to be better way - probably something to do
with reseting the plot region, but I can't recall how to do that now. If
there is, I'm sure someone will tell me what I overlooked.

Is this what you were looking for?

G

 
 Best,
 
 Georg
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

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


Re: [R] box-constrained

2008-03-05 Thread Ravi Varadhan
Hi,

Let me make the following points in response to your questions:

1.  Your call to optim() with L-BFGS-B as the method is correct.  Just
make sure that your function f is defined as negative log-likelihood,
since optim is by default a minimizer.  The other option is to define
log-likelihood as usual, but set control=list(fnscale=-1).

2.  You can add derivative (or gradient to be more precise) by defining that
function and then using the gr argument in optim.  Specifying exact
gradient almost always improves the convergence of the iterative schemes,
especially for ill-conditioned problems (flat region around the local
minima). So, if it is not too much trouble, and you are confident of your
differentiation skills, you should do that.  However, in most cases the
approximate finite-difference gradient used by optim() should be good
enough.

3.  Regardless of whether it is easy to compute the exact gradient or not,
it is generally a bad idea to maximize the likelihood that involves the
product of a large number of very small numbers.  It is almost always better
to maximize the log-likelihood.  Since the objective function is additive
rather than multiplicative, it has better numerical conditioning.

Ravi.



---

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: [EMAIL PROTECTED]

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

 





-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Gustave Lefou
Sent: Wednesday, March 05, 2008 1:34 PM
To: r-help@r-project.org
Subject: [R] box-constrained

Hello everybody,

I have a question about box-constrained optimization. I've done some
research and I found that optim could do that. Are there other ways in R ?

Is the following correct if I have a function f of two parameters belonging
for example to [0,1] and [0,Infinity] ?
optim(par=param, fn=f, method=L-BFGS-B, lower=c(0,0), upper=c(1,Inf))

My other question is whether it is possible to add the derivatives of my
function (like in nlm) and whether it is better to add them ?

If there is no need to add the derivatives, then I guess I could wish to
optimize the likelihood directly rather than to optimize the
log-likelihood... Indeed one aspect of the log-likelihood is to make the
derivatives tractable (in iid cases). Do you agree ?

Thank you !
Gustave

[[alternative HTML version deleted]]

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

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


Re: [R] N-way ANOVA

2008-03-05 Thread Prof Brian Ripley
On Wed, 5 Mar 2008, Paul Smith wrote:

 Dear All,

 Can R perform n-way ANOVA, i.e., with 3 or more factors?

Yes.  There are even examples on the help page!


 Thanks in advance,

 Paul

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


-- 
Brian D. Ripley,  [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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] inheritence in S4

2008-03-05 Thread Martin Morgan
[EMAIL PROTECTED] writes:

 Well well well...

You're partly misunderstanding...

 To summarize : let assume that A is a class (slot x) and C is a class
 containing A (A and slot y) - as(c,A) calls new(A). So new(A)
 HAS TO works, you can not decide to forbid empty object (unless you
 define setAs(C,A) ?)

new(A) has to work (return a valid object), but you can still create
arbitrary objects with a prototype (provided the prototype is
consistent with your validity function) and / or with initialize
methods (which do not have any named arguments, other than .Object,
without default values). Creating valid objects seems like a good
idea!

 - In addition, any test that you would like to set in initialize or
   validity have to first check is some field are empty (because
   if([EMAIL PROTECTED] 0)' will fail if x=numerical(0))

The validity and initialize methods have to be written in a way
consistent with valid objects -- if you want your object to contain
zero-length vectors, then the validity test has to accommodate that
(e.g., if (length([EMAIL PROTECTED])==0 || all([EMAIL PROTECTED]0)) {}). If 
your
initialize method is supposed to return a valid object, then you'd
better construct one! An initialize method on A that expected an
argument xValue whose default was 0 might be written as

setMethod(A, initialize, function(.Object, ..., xValue=0) {
callNextMethod(.Object, ..., x=xValue)
})

For this simple case you could create the object with a suitable
prototype and not have an initialize method at all

setClass(A, representation(x=numeric),
prototype=prototype(x=0))

 - When you call new(C), the neither new(A) nor intialize(A) are
   called. ()

new(A) is not called (why would it be?). If you have an initialize
method on A then it will be called. The problem you experienced
before was that your initialize method for A REQUIRED additional
arguments.

 So, all the security test you write in initialize A, you have to
 rewrite them on C ?

That is not correct; you do not have to duplicate code to construct an
object of class C.

 I start to undestand why not that much people use S4...

 Christophe

 [EMAIL PROTECTED] writes:

 Hi Martin, thanks for your answer

 But a couple
 of other quick points. I would have written

 setMethod(initialize, A,
  function(.Object, ..., xValue=numeric(0)){
  callNextMethod(.Object, ..., x=xValue)
  })


 I am not that much familiar with S3... In our way of writing this
 method, 'initialize' for 'A' will call the next method for A, which is
 initialise' for 'numeric', is that right ?

 from below,

 setClass(A, representation(x=numeric))

 The 'next' method refers to arguments in the signature of the generic,
 i.e., the 'next' method relevant to .Object. For A, the next method
 is the default initialize method, which will use all named arguments
 to fill the slots of .Object (and perhaps call validObject, among
 other things). You can see the code with

 showMethods(initialize, classes=ANY, includeDef=TRUE)

 S3 is not relevant here.

 And finally, the position of 'xValue' and 'yValue' means that
 the arugment has to be named, e.g., new(B, yValue=12). This seems a
 little awkward at first, but seems like a best practice

 I agree with you. But I do not like the use of ... , it lets us to
 make many mistake like in :

 print(3.5165,digitts=2)

 There is a typo in digitts, but since print works with ... , R does
 not detect the mistake.
 So I agree with the importance of naming argument, I always do it but
 without 'forcing' it

 If this were initialize, and you had provided an invalid named
 argument, the default method would have failed because there is no
 slot of that name.

 And finally, in Jim's thread I mention using a constructor. So in
 practice for a case like the above I would not define any initialize
 methods,

 Interesting. Will you still define a validity method or not even ?

 For a simple case like this there is no extra validity to check beyond
 that ensured by S4 (e.g., slots of the correct type). If there were
 constraints, e.g., only positive values, or length one vectors, then
 I would define a validity function.

 Martin


 Christophe



 Christophe Genolini [EMAIL PROTECTED] writes:

 Thanks Martin

 Well it works except that as seems to not like the initialize
 method : the following code (that is the same than yours with some
 initialize for A B and C) does not compile. It seems that as(c,A)
 does not work if we definie a initialize for A...

 --- 8 --
 setClass(A, representation(x=numeric))
 setMethod(initialize,A,function(.Object,value)[EMAIL PROTECTED] -
 value;return(.Object)})
 a - new(A,4)

 setClass(B, representation(y=numeric))
 setMethod(initialize,B,function(.Object,value)[EMAIL PROTECTED] -
 value;return(.Object)})
 b - new(B,5)

 setClass(C, contains=c(A, B))
 setMethod(initialize,C,function(.Object,valueA, valueB){
 [EMAIL PROTECTED] - valueA
 [EMAIL PROTECTED] - valueB
 

Re: [R] R_alloc with structures with flexible array members

2008-03-05 Thread Prof Brian Ripley
On Wed, 5 Mar 2008, Ramon Diaz-Uriarte wrote:

 Dear Jeff,

 Thanks for the suggestion. However, something is still not working.
 This is a simple example:

 ***  start C 
 #include R.h

 struct Sequence {
  int len;
  unsigned int state_count[];
 };


 int main(void) {

  struct Sequence *A;
  int n = 4;

  // First line segfaults. Second doesn't
A = (struct Sequence *) R_alloc(1,  sizeof(struct Sequence) + n *
 sizeof(unsigned int));
  //  A = malloc(sizeof(struct Sequence) + n * sizeof(unsigned int));

  return(0);
 }

 ***  end C  **


 I then do
 gcc -std=gnu99 -Wall -I/usr/share/R/include -I/usr/share/R/include
 -L/usr/lib/R/lib -lR ex7.c

 and the ./a.out segfaults when I use R_alloc (not with malloc).

You can't use R_alloc in a standalone program without initializing R, 
which has not been done here.

You said 'in a package', but this is not in a package.



 Best,

 R.

 On Wed, Mar 5, 2008 at 5:23 PM, Jeffrey Horner
 [EMAIL PROTECTED] wrote:
 Ramon Diaz-Uriarte wrote on 03/05/2008 04:25 AM:


 Dear All,
 
  In a package, I want to use some C code where I am using a structure
  (as the basic element of a linked list) with flexible array members.
  Basically, this is a structure where the last component is an
  incomplete array type  (e.g., Harbison  Steel, C, a reference
  manual, 5th ed., p. 159) such as:
 
  struct Sequence {
struct Sequence *next;
int len;
unsigned int state_count[];
  };
 
 
  To create one such sequence, I allocate storage (following Harbison
  and Steel) in a C program as follows:
 
  struct Sequence *A;
  int n = 4;
  A = malloc( sizeof(struct Sequence) + n * sizeof(unsigned int));
 
 
  If I understand correctly, however, it would be better to use R_alloc
  instead of malloc (memory automagically freed on exit and error;
  error-checking). But I do not know how to make the call to R_alloc
  here, since R_alloc allocates n units of size bytes each.
 
 
  I've tried, without success, the following two:
 
  int to_add_for_R_alloc =
  (int) ceil((float) sizeof(struct sequence) / sizeof(unsigned int));
 
A = (struct sequence *) R_alloc(to_add_for_R_alloc + n,
  sizeof(unsigned int));
 
  or even a brute force attempt as:
 
   A = (struct sequence *) R_alloc( 100,  sizeof(struct sequence));
 
 
  but both result in segmentation faults.
 
 
  Should I just keep using malloc (and free at end)?

  Hi Ramon,

  You should be able to use R_alloc without seg faults, so there's
  something wrong with your code somewhere. R_alloc multiplies its
  arguments together to come up with the total number of bytes to allocate
  then it allocates a raw vector and returns the data portion.

  So you can just treat R_alloc similarly to malloc by calling
  R_alloc(1,sizeof(struct Sequence) + n * sizeof(unsigned int)).

  Best,

  Jeff
  --
  http://biostat.mc.vanderbilt.edu/JeffreyHorner




 --
 Ramon Diaz-Uriarte
 Statistical Computing Team
 Structural Biology and Biocomputing Programme
 Spanish National Cancer Centre (CNIO)
 http://ligarto.org/rdiaz

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


-- 
Brian D. Ripley,  [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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] types of vectors / lists

2008-03-05 Thread a9804814
Hello,

I am an advanced user of R. Recently I found out that apparently I do 
not fully understand vectors and lists fully
Take this code snippet:


T = c(02.03.2008 12:23, 03.03.2008 05:54)
Times = strptime(T, %d.%m.%Y %H:%M)
Times # OK
class(Times)  # OK
is.list(Times)# sort of understand and not understand that
length(Times) # 9 ??? why is it length(Times[1]) ??

Times[1] # OK
Times[[1]] # Wrong


so Times is a vector-style thing of a non-atomic type.
What puzzles me is first that is.list returns true, and more importantly 
that the length-query returns the length of the first object of that 
apparent list and not how many elements are contained (i.e., I would 
have expected 2 as result - and I do not know what syntax to use in 
order to get the 2).

Moreover, if the whole thing is part of a data.frame:

DFTimes = as.data.frame(Times)
dim(DFTimes)
length(DFTimes$Times)  # OK, 2

then everything is as expected.

Could anyone please clearify why is.list returns true and why length in 
the first example returns 9 ? Is it that c() makes a list if the objects 
to be concatenated are not representable directly by atomic types ?

thanks,
Thomas

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


Re: [R] nls: different results if applied to normal or linearized data

2008-03-05 Thread Rolf Turner

On 6/03/2008, at 2:53 AM, Wolfgang Waser wrote:

 Dear all,

 I did a non-linear least square model fit

 y ~ a * x^b

 (a)  nls(y ~ a * x^b, start=list(a=1,b=1))

 to obtain the coefficients a  b.

 I did the same with the linearized formula, including a linear model

 log(y) ~ log(a) + b * log(x)

 (b)  nls(log10(y) ~ log10(a) + b*log10(x), start=list(a=1,b=1))
 (c)  lm(log10(y) ~ log10(x))

 I expected coefficient b to be identical for all three cases.  
 Hoever, using my
 dataset, coefficient b was:
 (a) 0.912
 (b) 0.9794
 (c) 0.9794

 Coefficient a also varied between option (a) and (b), 107.2 and 94.7,
 respectively.

 Is this supposed to happen? Which is the correct coefficient b?

The two approaches assume two different models.

Model (1) is y = a*x^b + E (where different errors are independent  
and identically
--- usually normally --- distributed).

Model (2) is y = a*(x^b)*E (and you are usually tacitly assuming  
that ln E is normally distributed).

The point estimates of a and b will consequently be different ---  
although usually not hugely
different.  Their distributional properties will be substantially  
different.

cheers,

Rolf Turner

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

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


[R] collapsing list to data.frame

2008-03-05 Thread a9804814
Hello,


Given a list with all elements having identical layout, e.g.:


l = NULL
l[[1]] = list(4, hello)
l[[2]] = list(7, world)
l[[3]] = list(9,   )


is there an easy way to collapse this list into a data.frame with each 
row being the elements of the list ?
I.e. in this case I want to convert the list into a data.frame with 3 
rows and 2 columns, where column 1 holds the integer values, and column 
2 the character values.

I can get it done by looping over all elements and rbind them together 
to the final result, but that is quite slow (for large sets) and ugly, 
so I was wondering if there's an easy syntax.

thanks a lot in advance,
Thomas

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


Re: [R] N-way ANOVA

2008-03-05 Thread Paul Smith
On Wed, Mar 5, 2008 at 7:41 PM, Prof Brian Ripley [EMAIL PROTECTED] wrote:
   Can R perform n-way ANOVA, i.e., with 3 or more factors?

  Yes.  There are even examples on the help page!

Thanks, Prof. Ripley. I will have a look at it.

Paul

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


[R] differentiating a numeric vector

2008-03-05 Thread Levi Waldron
What functions exist for differentiating a numeric vector (in my case
spectral data)?  That is, experimental data without an analytical
function.  ie,

 x - seq(1,10,0.1)
 y=x^3+rnorm(length(x),sd=0.01) #although the real function would be 
 nothing simple like x^3...
 derivy - 

I know I could just use diff(y) but it would be nice to estimate
derivatives at the endpoints, calculate higher-order derivatives, and
maybe have some smoothing options ie by differentiating a spline or
something like that.

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


Re: [R] legend for several graphics

2008-03-05 Thread hadley wickham
On Wed, Mar 5, 2008 at 8:28 AM, Georg Otto [EMAIL PROTECTED] wrote:
 Hi,

  I am trying to generate a figure of 9 plots that are contained in one
  device by using

  par(mfrow = c(3,3,))

  I would like to have 1 common legend for all 9 plots somewhere outside
  of the plotting area (as opposed to one legend inside each of the 9
  plots, which the function legend() seems to generate by default).

If you provide more detail about your problem (what are the 9 plots?)
I'm sure we can suggest other solutions using lattice or ggplot that
will substantially simplify your code, as well as automatically
drawing the legend.

Hadley


-- 
http://had.co.nz/

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


Re: [R] Asking, are simple effects different from 0

2008-03-05 Thread jebyrnes

Indeed, but are not each of the cell means also evaluations of the effect of
one factor at the specific level of another factor?  Is this an issue of
Tomato, tomahto.

I guess my question is, if I want to know if each of those is different from
0, then should I use the 48df from the full model, or the 9 for each cell?


Chuck Cleland wrote:
 
That does not corresponds to what I think of as the simple effects. 
 That specifies the six cell means, but it does not *compare* any cell 
 means.  I think of a simple effect as the effect of one factor at a 
 specific level of some other factor.
 
 summary(glht(fm, linfct = cm2), test = adjusted(type=none)) 
 
 Correct? What is the df on those t-tests then?  Is it 48?
 
Yes, df = 48 for each contrast.
 
 Interestingly, I find this produces results no different than
 
 fm2-lm(breaks ~ tension:wool+0, data=warpbreaks) 
 summary(fm2)
 
Yes, but those are not what I would call the simple effects.  Those 
 are essentially one-sample t-tests for each of the 6 cell means.
 
 Also, here, it would seem each t-test was done with the full 48df.  Hrm.
 
The df are based on the whole model, not the 9 observations in one
 cell.
 
 
 

-- 
View this message in context: 
http://www.nabble.com/Asking%2C-are-simple-effects-different-from-0-tp15835552p15857771.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] collapsing list to data.frame

2008-03-05 Thread Erik Iverson
Perhaps

data.frame(do.call(rbind, l))

?

- Erik Iverson

[EMAIL PROTECTED] wrote:
 Hello,
 
 
 Given a list with all elements having identical layout, e.g.:
 
 
 l = NULL
 l[[1]] = list(4, hello)
 l[[2]] = list(7, world)
 l[[3]] = list(9,   )
 
 
 is there an easy way to collapse this list into a data.frame with each 
 row being the elements of the list ?
 I.e. in this case I want to convert the list into a data.frame with 3 
 rows and 2 columns, where column 1 holds the integer values, and column 
 2 the character values.
 
 I can get it done by looping over all elements and rbind them together 
 to the final result, but that is quite slow (for large sets) and ugly, 
 so I was wondering if there's an easy syntax.
 
 thanks a lot in advance,
 Thomas
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


[R] read.table

2008-03-05 Thread Georg Ehret
Dear R community,   I encounter a problem reading data into a dataframe. See
attachment for the input. I use:
data-read.table(test,fill=T,row.names=1)

When you look at data you can see that some lines of my input were
broken... I can avoid this problem by sorting the lines by length... Do I
have it wrong or is there a bug?

Thanks and wishing you a great remainder of the day,
Georg.
**
Georg Ehret
JHU
Baltimore - US
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] collapsing list to data.frame

2008-03-05 Thread Dimitris Rizopoulos
try this:

l - vector(list, 3)
l[[1]] - list(4, hello)
l[[2]] - list(7, world)
l[[3]] - list(9,   )

lis - lapply(l, names-, value = c(V1, V2))
do.call(rbind, lapply(lis, data.frame, stringsAsFactors = FALSE))

I hope it helps.

Best,
Dimitris


Dimitris Rizopoulos
Biostatistical Centre
School of Public Health
Catholic University of Leuven

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


Quoting [EMAIL PROTECTED]:

 Hello,


 Given a list with all elements having identical layout, e.g.:


 l = NULL
 l[[1]] = list(4, hello)
 l[[2]] = list(7, world)
 l[[3]] = list(9,   )


 is there an easy way to collapse this list into a data.frame with each
 row being the elements of the list ?
 I.e. in this case I want to convert the list into a data.frame with 3
 rows and 2 columns, where column 1 holds the integer values, and column
 2 the character values.

 I can get it done by looping over all elements and rbind them together
 to the final result, but that is quite slow (for large sets) and ugly,
 so I was wondering if there's an easy syntax.

 thanks a lot in advance,
 Thomas

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





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

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


[R] Axes in polymars

2008-03-05 Thread Shewcraft, Ryan
Hi All,

I can't quite figure out how to change the parameters of the x and y
axes when I plot a polymars object.  I want to add a scatterplot of the
data points, but the polymars plot seems to automatically set the
parameters to fit the polymars line.  I've tried using plot(poly,1,fig=
c(x1,x2,y1,y2)) but have no luck.  Any thoughts?

Thanks,
Ryan 
 


Please open the following attachment for important information regarding this 
e-mail communication.
This email communication is privileged, confidential or otherwise protected
by disclosure and is intended only for the individuals or entities named
above and any others who have been specifically authorized to receive it. 
Any unauthorized dissemination, copying or use of the contents of this
email is strictly prohibited and may be in violation of law.  If you are
not the intended recipient, please do not read, copy, use or disclose to
others the contents of this communication.  Please notify the sender that
you have received this e-mail in error by replying to this e-mail or by
phoning +1-212-583-5000 Monday to Friday between 7:30 am and 7:30 pm (EST).
At any other time, please call +1-212-583-5221.  Please then delete the
e-mail from your system and any copies of it.  Nothing contained in this
disclaimer shall be construed in any way to grant permission to transmit
confidential information via this firm's e-mail system or as a waiver of
any confidentiality or privilege.
 
This communication may contain highly confidential information regarding
Blackstone's investments, strategy and organization.  Your acceptance of
this communication from Blackstone constitutes your agreement to (i) keep
confidential all the confidential information contained in this
communication, as well as any information derived by you from the
confidential information contained in this communication (collectively,
Confidential Information) and not disclose any such Confidential
Information to any other person, (ii) not use any of the Confidential
Information for any purpose other than pursuant to your relationship with
Blackstone, (iii) not use the Confidential Information for purposes of
trading any security, including, without limitation, securities of
Blackstone or its portfolio companies, (iv) not copy any Confidential
Information without Blackstone's prior consent and (v) promptly return any
Confidential Information contained in this communication to Blackstone upon
Blackstone's request.
 
This communication is for informational purposes only and does not
constitute an offer to sell or a solicitation of an offer to purchase any
interest in any investment vehicles managed by any business unit of The
Blackstone Group.  Blackstone does not accept any responsibility or
liability arising from the use of this communication.  No representation is
being made that the information presented is accurate, current or complete,
and such information is at all times subject to change without notice. 
Opinions expressed may differ or be contrary to the opinions and
recommendations of a Blackstone business unit.  
 
Blackstone does not provide legal, accounting or tax advice.  Any statement
regarding legal, accounting or tax matters was written in connection with
the explanation of the matters described herein and was not intended or
written to be relied upon by any person as definitive advice.  Any
discussion of U.S. tax matters contained within this communication is not
intended to be used and cannot be used for the purpose of avoiding 
penalties that may be imposed under applicable Federal, state or local tax
law or recommending to another party any transaction or matter addressed
herein.  Each person should seek advice based on its particular
circumstances from independent legal, accounting and tax advisors regarding
the matters discussed in this e-mail.__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] types of vectors / lists

2008-03-05 Thread Erik Iverson
Hello -

[EMAIL PROTECTED] wrote:
 Hello,
 
 I am an advanced user of R. Recently I found out that apparently I do 
 not fully understand vectors and lists fully
 Take this code snippet:
 
 
 T = c(02.03.2008 12:23, 03.03.2008 05:54)
 Times = strptime(T, %d.%m.%Y %H:%M)
 Times # OK
 class(Times)  # OK
 is.list(Times)# sort of understand and not understand that

?POSIXlt says

Class 'POSIXlt' is a named list of
  vectors representing

  'sec' 0-61: seconds

  'min' 0-59: minutes

  'hour' 0-23: hours

  'mday' 1-31: day of the month

  'mon' 0-11: months after the first of the year.

  'year' Years since 1900.

  'wday' 0-6 day of the week, starting on Sunday.

  'yday' 0-365: day of the year.

  'isdst' Daylight savings time flag.  Positive if in force, zero if
   not, negative if unknown.


 length(Times) # 9 ??? why is it length(Times[1]) ??

Because Times is a list with 9 elements, try typing names(Times) to see 
what those elements are called.  That will help you understand.

 
 Times[1] # OK
 Times[[1]] # Wrong

Whether or not this is 'wrong' has been I believe debated before on the 
list.  See this message and the rest of the thread associated with it.

http://tolstoy.newcastle.edu.au/R/help/06/07/31508.html

- Erik Iverson

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


  1   2   >