[R] Hand-crafting an .RData file

2009-11-08 Thread Adam D. I. Kramer

Hello,

I frequently have to export a large quantity of data from some
source (for example, a database, or a hand-written perl script) and then
read it into R.  This occasionally takes a lot of time; I'm usually using
read.table("filename",comment.char="",quote="") to read the data once it is
written to disk.

However, I *know* that the program that generates the data is more
or less just calling printf in a for loop to create the csv or tab-delimited
file, writing, then having R parse it, which is pretty inefficient. 
Instead, I am interested in figuring out how to write the data in .RData

format so that I can load() it instead of read.table() it.

Trolling the internet, however, has not suggested anything about the
specification for an .RData file. Could somebody link me to a specification
or some information that would instruct me on how to construct a .RData
file (either compressed or uncompressed)?

Also, I am open to other suggestions of how to get load()-like
efficiency in some other way.

Many thanks,
Adam D. I. Kramer

__
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] How to change color the default in levelplot() ?

2009-11-08 Thread BabyL BabyK
Dear R communities

May I seek your advices on how to change color the default in levelplot(), e.g. 
from the default of pink and light blue, to e.g. red and green ? 
The levelplot function has 1 of the arguments being panel (which is actually 
panel.levelplot), but I am not sure where the commands to alter the color. 

For example, I type:
p1<-levelplot(my.mat,colorkey=FALSE),
how could I include col.regions of panel.levelplot? 
And whether it is right to use the col.regions? 

I am using R 2.8.1 in ubuntu.

Many thanks and have a good day! 



  
[[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] Forward selection peculiarity

2009-11-08 Thread ricefan81

Hi all,

I'm using the mle.stepwise function in the wle package to run forward
selection. My statistical background tells me that the first variable
selected should be the one with the highest correlation with the response,
however that's not the case. The two highest correlations with the response
are "similar" (0.86 and 0.89) however the forward selection selects the
variable with 0.86 correlation first and then in the next step chooses the
variable with correlation 0.89.

Any thoughts on this? Is is this some type of bug in the function or
package?

Thanks
-- 
View this message in context: 
http://old.nabble.com/Forward-selection-peculiarity-tp26260670p26260670.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] Outputing multilple subsets

2009-11-08 Thread Johann Hibschman

On Nov 8, 2009, at 7:23 PM, rusers.sh wrote:


for (i in num)  {
  c_num<-c[c$b==num,]
  write.csv(c_num,file="c:/c_num.csv")
}

Warning messages:
1: In c$b == num :
 longer object length is not a multiple of shorter object length


This is because you're comparing column b to the entire vector of  
numbers (num), not the current number in the iteration (i).  The first  
line of the loop should be "c_num<-c[c$b==i,]".


From a style point of view, I'd use "n" as my variable, since "i" is  
too commonly used as an integer index.


Also, you will be overwriting the same file, called "c_num.csv", on  
each iteration.


You should try something more like:

for (n in num) {
  c.n <- c[c$b==n,]
  write.csv(c.n, file=paste("c:/c_", n, ".csv", sep="")
}

I hope that helps.

Cheers,
Johann Hibschman

__
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] Bug in gplots::heatmap.2 symbreaks arg default (?)

2009-11-08 Thread Steve Lianoglou

Hi all,

I think the code to calculate the default value of the symbreaks  
argument in the gplots::heatmap.2 function is wrong.


The documentation says symbreaks defaults to true if any negative  
values are detected in the data passed in.


The relevant code in the parameter list of this function definition  
(gplots 2.7.3) is this:


symbreaks = min(x < 0, na.rm = TRUE) || scale != "none"

When I'm pretty sure it should be:

symbreaks = min(x, na.rm = TRUE) < 0 || scale != "none"

Likewise for the symkey parameter.

Right?

Thanks,
-steve

--
Steve Lianoglou
Graduate Student: Computational Systems Biology
  |  Memorial Sloan-Kettering Cancer Center
  |  Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact

__
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] reference on contr.helmert and typo on its help page.

2009-11-08 Thread John Fox
Dear Peng,

I'm tempted to try to get an entry in the fortunes package but will instead
try to answer your questions directly:

> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
> Behalf Of Peng Yu
> Sent: November-08-09 7:41 PM
> To: r-h...@stat.math.ethz.ch
> Subject: Re: [R] reference on contr.helmert and typo on its help page.
> 
> Dear John,
> 
> I did read Section 9.1.2 and various other textbooks before posting my
> questions. But each reference uses slightly different notations and
> terminology. I get confused and would like a description that
> summaries everything so that I don't have to refer to many different
> resources. May I ask a few questions on the section in your textbook?
> 
> Which variable in Section 9.1.2 is "a matrix of contrasts" mentioned
> in the help page of 'contr.helmert'? Which matrix of contrast in R
> corresponds to dummy regression? With different R formula, e.g. y ~ x
> vs. y ~ x -1, $X_F$ (mentioned on page 189) is different and hence
> $\beta_F$ (mentioned in eq. 9.3) is be different. So my understanding
> is that the matrix of contrast should depend on the formula. But it is
> not according to the help page of "contr.helmert".

If the model is simply y ~ A, for the factor A, then cbind(1, contrasts(A))
is what I call X_B, the row-basis of the model matrix. As I explain in the
section that you read, the level means are mu = X_B beta, and thus beta =
X_B^-1 mu = 0 are the hypotheses tested by the contrasts. Moreover, if, as
in Helmert contrasts, the columns of X_B are orthogonal, then so are the
rows of X_B^-1, and the latter are simply rescalings of the former. That
allows one conveniently to code the hypotheses directly in X_B; all this is
also explained in that section of my book, and is essentially what Peter D.
told you. In R, contr.treatment and contr.SAS provide dummy-variable (0/1)
coding of regressors, differing only in the selection of the reference
level.

John

> 
> Regards,
> Peng
> 
> On Sun, Nov 8, 2009 at 6:17 PM, John Fox  wrote:
> > Dear Peng Yu,
> >
> > Perhaps you're referring to my text, Applied Linear Regression Analysis
and
> > Generalized Linear Models, since I seem to recall that you sent me a
number
> > of questions about it. See Section 9.1.2 on linear contrasts for the
answer
> > to your question.
> >
> > I hope this helps,
> >  John
> >
> > 
> > John Fox
> > Senator William McMaster
> >  Professor of Social Statistics
> > Department of Sociology
> > McMaster University
> > Hamilton, Ontario, Canada
> > web: socserv.mcmaster.ca/jfox
> >
> >
> >> -Original Message-
> >> From: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org]
> > On
> >> Behalf Of Peng Yu
> >> Sent: November-08-09 4:52 PM
> >> To: r-h...@stat.math.ethz.ch
> >> Subject: Re: [R] reference on contr.helmert and typo on its help page.
> >>
> >> On Sun, Nov 8, 2009 at 11:28 AM, Peter Dalgaard
> >>  wrote:
> >> > Gabor Grothendieck wrote:
> >> >>
> >> >> On Sun, Nov 8, 2009 at 11:59 AM, Peng Yu 
wrote:
> >> >>>
> >> >>> On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch

> >> >>> wrote:
> >> 
> >>  On 08/11/2009 11:03 AM, Peng Yu wrote:
> >> >
> >> > I'm wondering which textbook discussed the various contrast
matrices
> >> > mentioned in the help page of 'contr.helmert'. Could somebody let
me
> >> > know?
> >> 
> >>  Doesn't the reference on that page discuss them?
> >> >>>
> >> >>> It does explain what the functions are. But I need a more basic and
> >> >>> complete reference. For example, I want to understand what 'Helmert
> >> >>> parametrization' (on page 33 of 'Statistical Models in S') is.
> >> >>>
> >> >>
> >> >> Just google for: Helmert contrasts
> >> >
> >> > Or,
> >> >
> >> >> contr.helmert(5)
> >> >  [,1] [,2] [,3] [,4]
> >> > 1   -1   -1   -1   -1
> >> > 2    1   -1   -1   -1
> >> > 3    0    2   -1   -1
> >> > 4    0    0    3   -1
> >> > 5    0    0    0    4
> >> >
> >> >> MASS::fractions(MASS::ginv(contr.helmert(5)))
> >> >     [,1]  [,2]  [,3]  [,4]  [,5]
> >> > [1,]  -1/2   1/2     0     0     0
> >> > [2,]  -1/6  -1/6   1/3     0     0
> >> > [3,] -1/12 -1/12 -1/12   1/4     0
> >> > [4,] -1/20 -1/20 -1/20 -1/20   1/5
> >> >
> >> > and apply brains.
> >> >
> >> > I.e., except for a slightly odd multiplier, the parameters represent
the
> >> >  difference between each level and the average of the preceding
levels.
> >>
> >> I realized that my questions are what a contrast matrix is and how it
> >> is related to hypothesis testing. For a give hypothesis, how to get
> >> the corresponding contrast matrix in a systematical way? There are
> >> some online materials, but they are all diffused. I have also read the
> >> book Applied Linear Regression Models, which doesn't give a complete
> >> descriptions on all the aspects of contrast and contrast matrix. But I
> >> would want a textbook that gives a complete description, so that I

Re: [R] Models for Discrete Choice in R

2009-11-08 Thread Frank E Harrell Jr

Emmanuel Charpentier wrote:

Le dimanche 08 novembre 2009 à 17:07 -0200, Iuri Gavronski a écrit :

Hi,

I would like to fit Logit models for ordered data, such as those
suggested by Greene (2003), p. 736.

Does anyone suggests any package in R for that?


look up the polr function in package MASS (and read the relevant pages
in V&R4 and some quoted references...) or the slightly more
sophisticated (larger range of models) lrm function in F. Harrell's
Design (now rms) packge (but be aware that Design is a huge beast witch
carries its own "computing universe", based on (strong) Harrell's view
of what a regression analysis should be : reading his book is, IMHO,
necessary to understand his choices and agree (or disgree) with them).

If you have a multilevel model (a. k. a. one "random effect" grouping),
the "repolr" packge aims at that, but I've been unable to use it
recently (numerical exceptions).


By the way, my dependent variable is ordinal and my independent
variables are ratio/intervalar.


Numeric ? Then maybe some recoding/transformation is in order ... in
which case Design/rms might or might not be useful.


I'm not clear on what recoding or transformation is needed for an 
ordinal dependent variable and ratio/interval independent variables, nor 
why rms/Design would not be useful.


Frank



HTH,

Emmanuel Charpentier



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

__
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] Outputing multilple subsets

2009-11-08 Thread Ista Zahn
Have you considered using split?

-Ista

On Sun, Nov 8, 2009 at 7:23 PM, rusers.sh  wrote:
> Hi Rusers,
>  I hope to divide the original dataset into several subsets and output
> these multilple datasets. But errors appeared in my loops. See example.
> ##
> a<-c(1:10)
> b<-c(rep(1,3),rep(2,3),rep(3,4))
> c<-data.frame(a,b)  #c is the example data
> num<-c(unique(b))
> # I hope to get the subsets c_1.csv,c_2.csv and c_3.csv
> #Errors
> for (i in num)  {
>   c_num<-c[c$b==num,]
>   write.csv(c_num,file="c:/c_num.csv")
> }
>
> Warning messages:
> 1: In c$b == num :
>  longer object length is not a multiple of shorter object length
> 2: In c$b == num :
>  longer object length is not a multiple of shorter object length
> 3: In c$b == num :
>  longer object length is not a multiple of shorter object length
>  I think the problem should be file="c:/c_num.csv", anybody has ever met
> this problem?
>  Thanks very much.
> -
> Jane Chang
> Queen's
>
>        [[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.
>



-- 
Ista Zahn
Graduate student
University of Rochester
Department of Clinical and Social Psychology
http://yourpsyche.org

__
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] Simple 2-Way Anova issue in R

2009-11-08 Thread jim holtman
> x
  Year Depth Biomass1 Biomass2
1 199910 14.3 14.7
2 199915 14.7 15.6
> require(reshape)
> melt(x, id=c('Year','Depth'))
  Year Depth variable value
1 199910 Biomass1  14.3
2 199915 Biomass1  14.7
3 199910 Biomass2  14.7
4 199915 Biomass2  15.6


On Sun, Nov 8, 2009 at 7:38 PM, znd  wrote:
>
> Ah, I believe I constructed my *.csv wrong in that I only had 1 observation
> within groups whereas I needed at least 2.
>
> Originally I had:
>
> Year     Depth     Biomass1     Biomass2
> 1999     10         14.3           14.7
> 1999     15         14.7           15.6
> etc.
>
> but I switched this to:
>
> Year     Depth     Biomass
> 1999     10         14.3
> 1999     10         14.7
> 1999     15         15.1
> 1999     15         15.6
>
> I believe this may be the appropriate way to collapse my biomass data in
> order to perform a 2-way anova.
>
> If not feel free comment and thanks for the help.
> --
> View this message in context: 
> http://old.nabble.com/Simple-2-Way-Anova-issue-in-R-tp26258684p26259649.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.
>



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

What is the problem that 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.


Re: [R] reference on contr.helmert and typo on its help page.

2009-11-08 Thread Peng Yu
Dear John,

I did read Section 9.1.2 and various other textbooks before posting my
questions. But each reference uses slightly different notations and
terminology. I get confused and would like a description that
summaries everything so that I don't have to refer to many different
resources. May I ask a few questions on the section in your textbook?

Which variable in Section 9.1.2 is "a matrix of contrasts" mentioned
in the help page of 'contr.helmert'? Which matrix of contrast in R
corresponds to dummy regression? With different R formula, e.g. y ~ x
vs. y ~ x -1, $X_F$ (mentioned on page 189) is different and hence
$\beta_F$ (mentioned in eq. 9.3) is be different. So my understanding
is that the matrix of contrast should depend on the formula. But it is
not according to the help page of "contr.helmert".

Regards,
Peng

On Sun, Nov 8, 2009 at 6:17 PM, John Fox  wrote:
> Dear Peng Yu,
>
> Perhaps you're referring to my text, Applied Linear Regression Analysis and
> Generalized Linear Models, since I seem to recall that you sent me a number
> of questions about it. See Section 9.1.2 on linear contrasts for the answer
> to your question.
>
> I hope this helps,
>  John
>
> 
> John Fox
> Senator William McMaster
>  Professor of Social Statistics
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> web: socserv.mcmaster.ca/jfox
>
>
>> -Original Message-
>> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
> On
>> Behalf Of Peng Yu
>> Sent: November-08-09 4:52 PM
>> To: r-h...@stat.math.ethz.ch
>> Subject: Re: [R] reference on contr.helmert and typo on its help page.
>>
>> On Sun, Nov 8, 2009 at 11:28 AM, Peter Dalgaard
>>  wrote:
>> > Gabor Grothendieck wrote:
>> >>
>> >> On Sun, Nov 8, 2009 at 11:59 AM, Peng Yu  wrote:
>> >>>
>> >>> On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch 
>> >>> wrote:
>> 
>>  On 08/11/2009 11:03 AM, Peng Yu wrote:
>> >
>> > I'm wondering which textbook discussed the various contrast matrices
>> > mentioned in the help page of 'contr.helmert'. Could somebody let me
>> > know?
>> 
>>  Doesn't the reference on that page discuss them?
>> >>>
>> >>> It does explain what the functions are. But I need a more basic and
>> >>> complete reference. For example, I want to understand what 'Helmert
>> >>> parametrization' (on page 33 of 'Statistical Models in S') is.
>> >>>
>> >>
>> >> Just google for: Helmert contrasts
>> >
>> > Or,
>> >
>> >> contr.helmert(5)
>> >  [,1] [,2] [,3] [,4]
>> > 1   -1   -1   -1   -1
>> > 2    1   -1   -1   -1
>> > 3    0    2   -1   -1
>> > 4    0    0    3   -1
>> > 5    0    0    0    4
>> >
>> >> MASS::fractions(MASS::ginv(contr.helmert(5)))
>> >     [,1]  [,2]  [,3]  [,4]  [,5]
>> > [1,]  -1/2   1/2     0     0     0
>> > [2,]  -1/6  -1/6   1/3     0     0
>> > [3,] -1/12 -1/12 -1/12   1/4     0
>> > [4,] -1/20 -1/20 -1/20 -1/20   1/5
>> >
>> > and apply brains.
>> >
>> > I.e., except for a slightly odd multiplier, the parameters represent the
>> >  difference between each level and the average of the preceding levels.
>>
>> I realized that my questions are what a contrast matrix is and how it
>> is related to hypothesis testing. For a give hypothesis, how to get
>> the corresponding contrast matrix in a systematical way? There are
>> some online materials, but they are all diffused. I have also read the
>> book Applied Linear Regression Models, which doesn't give a complete
>> descriptions on all the aspects of contrast and contrast matrix. But I
>> would want a textbook that gives a complete description, so that I
>> don't have to look around for other materials.
>>
>> __
>> 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] Simple 2-Way Anova issue in R

2009-11-08 Thread znd

Ah, I believe I constructed my *.csv wrong in that I only had 1 observation
within groups whereas I needed at least 2.

Originally I had:

Year Depth Biomass1 Biomass2
1999 10 14.3   14.7
1999 15 14.7   15.6
etc.

but I switched this to:

Year Depth Biomass
1999 10 14.3
1999 10 14.7
1999 15 15.1
1999 15 15.6

I believe this may be the appropriate way to collapse my biomass data in
order to perform a 2-way anova.

If not feel free comment and thanks for the help.
-- 
View this message in context: 
http://old.nabble.com/Simple-2-Way-Anova-issue-in-R-tp26258684p26259649.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] Outputing multilple subsets

2009-11-08 Thread rusers.sh
Hi Rusers,
  I hope to divide the original dataset into several subsets and output
these multilple datasets. But errors appeared in my loops. See example.
##
a<-c(1:10)
b<-c(rep(1,3),rep(2,3),rep(3,4))
c<-data.frame(a,b)  #c is the example data
num<-c(unique(b))
# I hope to get the subsets c_1.csv,c_2.csv and c_3.csv
#Errors
for (i in num)  {
   c_num<-c[c$b==num,]
   write.csv(c_num,file="c:/c_num.csv")
}

Warning messages:
1: In c$b == num :
  longer object length is not a multiple of shorter object length
2: In c$b == num :
  longer object length is not a multiple of shorter object length
3: In c$b == num :
  longer object length is not a multiple of shorter object length
  I think the problem should be file="c:/c_num.csv", anybody has ever met
this problem?
  Thanks very much.
-
Jane Chang
Queen's

[[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] reference on contr.helmert and typo on its help page.

2009-11-08 Thread John Fox
Dear Peng Yu,

Perhaps you're referring to my text, Applied Linear Regression Analysis and
Generalized Linear Models, since I seem to recall that you sent me a number
of questions about it. See Section 9.1.2 on linear contrasts for the answer
to your question.

I hope this helps,
 John


John Fox
Senator William McMaster 
  Professor of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
> Behalf Of Peng Yu
> Sent: November-08-09 4:52 PM
> To: r-h...@stat.math.ethz.ch
> Subject: Re: [R] reference on contr.helmert and typo on its help page.
> 
> On Sun, Nov 8, 2009 at 11:28 AM, Peter Dalgaard
>  wrote:
> > Gabor Grothendieck wrote:
> >>
> >> On Sun, Nov 8, 2009 at 11:59 AM, Peng Yu  wrote:
> >>>
> >>> On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch 
> >>> wrote:
> 
>  On 08/11/2009 11:03 AM, Peng Yu wrote:
> >
> > I'm wondering which textbook discussed the various contrast matrices
> > mentioned in the help page of 'contr.helmert'. Could somebody let me
> > know?
> 
>  Doesn't the reference on that page discuss them?
> >>>
> >>> It does explain what the functions are. But I need a more basic and
> >>> complete reference. For example, I want to understand what 'Helmert
> >>> parametrization' (on page 33 of 'Statistical Models in S') is.
> >>>
> >>
> >> Just google for: Helmert contrasts
> >
> > Or,
> >
> >> contr.helmert(5)
> >  [,1] [,2] [,3] [,4]
> > 1   -1   -1   -1   -1
> > 2    1   -1   -1   -1
> > 3    0    2   -1   -1
> > 4    0    0    3   -1
> > 5    0    0    0    4
> >
> >> MASS::fractions(MASS::ginv(contr.helmert(5)))
> >     [,1]  [,2]  [,3]  [,4]  [,5]
> > [1,]  -1/2   1/2     0     0     0
> > [2,]  -1/6  -1/6   1/3     0     0
> > [3,] -1/12 -1/12 -1/12   1/4     0
> > [4,] -1/20 -1/20 -1/20 -1/20   1/5
> >
> > and apply brains.
> >
> > I.e., except for a slightly odd multiplier, the parameters represent the
> >  difference between each level and the average of the preceding levels.
> 
> I realized that my questions are what a contrast matrix is and how it
> is related to hypothesis testing. For a give hypothesis, how to get
> the corresponding contrast matrix in a systematical way? There are
> some online materials, but they are all diffused. I have also read the
> book Applied Linear Regression Models, which doesn't give a complete
> descriptions on all the aspects of contrast and contrast matrix. But I
> would want a textbook that gives a complete description, so that I
> don't have to look around for other materials.
> 
> __
> 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] Express "e" in R?

2009-11-08 Thread David Winsemius

exp(1)
#[1] 2.718282

On Nov 8, 2009, at 5:03 PM, Crab wrote:



How do you express "e" the base of the natural log in R.
--



David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] Express "e" in R?

2009-11-08 Thread Duncan Murdoch

On 08/11/2009 5:03 PM, Crab wrote:
How do you express "e" the base of the natural log in R. 


e <- exp(1)

__
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] Express "e" in R?

2009-11-08 Thread Crab

How do you express "e" the base of the natural log in R. 
-- 
View this message in context: 
http://old.nabble.com/Express-%22e%22-in-R--tp26258503p26258503.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] Normal distribution

2009-11-08 Thread mandalic05

Normal distribution check within R can be done with functions available in
nortest package. This package consists of several normality tests. In order
to install package type install.packages('nortest'). Afterwards, you should
consider running ks.test() only if mu and sigma parameters are known (these
stand for population arithmetic mean and variance) - and that's only
applicable if your data is gathered from the population. Therefor I
recommend using lillie.test() function, which is Lilliefors' modification of
Kolmogorov-Smirnov statistic. It's applicable both for data gathered from a
sample, but can also be applied to population data. You can run ks.test(x,
pnorm), but don't worry if you get several ties - these occur due to
rounding of values, or if your data come from descrete probability
function...

You can also try shapiro.test() function if your sample counts less then 50
responses (Shapiro-Wilks' test for small samples), or ad.test() for
Anderson-Darling normality test. You should revise these statistical
procedures in official literature, but there's also a lot of info on
wikipedia about stated statistical techniques.

If absolute value of skewness is larger than 1.96 * standard error of
skewness, your distribution significantly differs from normal. Also stands
for kurtosis. Value 1.96 implies p-value lower than .05, and 2.58 lower than
.01
Skewness function is called with skew(), and kurtosis with kurtosi()
function.

Standard error of skewness is calculated from formula se.sk <-
sqrt(6/length(x)) and standard error of kurtosis from formula se.ku <-
sqrt(24/length(x))

If mean is not located in the middle of the range, this can also indicate a
violation of normality.

I strongly recommend reading official help for nortest package, and
consulting an official statistical literature!

P.S.

shapiro.test() is located in stats package, but running it on a large sample
(N >> 50) is not quite applicable, hence use lillie.test() for those
purposes, or ks.test(x, pnorm) - where x argument is variable which
normality you're about to check, and pnorm stands for normal distribution
function.



Noela Sánchez-2 wrote:
> 
> Hi,
> 
> I am dealing with how to check in R if some data that I have belongs to a
> normal distribution or not. I am not interested in obtaining the
> theoreticall frequencies. I am only interested in determing if (by means
> of
> a test as Kolmogorov, or whatever), if my data are normal or not.
> 
> But I have tried with ks.test() and I have not got it.
> 
> 
> -- 
> Noela
> Grupo de Recursos Marinos y Pesquerías
> Universidad de A Coruña
> 
>   [[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.
> 
> 

-- 
View this message in context: 
http://old.nabble.com/Normal-distribution-tp25702570p26259120.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] Turn dates into age

2009-11-08 Thread Marc Schwartz

> Sys.Date()
[1] "2009-11-08"

> as.Date("05/29/1971", format = "%m/%d/%Y")
[1] "1971-05-29"


> as.numeric((Sys.Date() - as.Date("05/29/1971", format = "%m/%d/ 
%Y")) / 365.25)

[1] 38.44764

or perhaps more clearly:

EndDate <- Sys.Date()
StartDate <- as.Date("05/29/1971", format = "%m/%d/%Y")

> as.numeric((EndDate - StartDate) / 365.25)
[1] 38.44764



We coerce to numeric here, to return a standard numeric value, rather  
than the result being of class difftime with an attribute of 'days'  
for units:


> str((Sys.Date() - as.Date("05/29/1971", format = "%m/%d/%Y")) /  
365.25)

Class 'difftime'  atomic [1:1] 38.4
  ..- attr(*, "units")= chr "days"



HTH,

Marc Schwartz


On Nov 8, 2009, at 5:22 PM, Jim Burke wrote:


To clarify.

Lets turn a date into an age. Given 05/29/1971 in mm/dd/
format. What is the year difference between then and today?
This would be the "age" requested that starts 05/29/1971 as
one.

Thanks,
Jim




David Winsemius wrote:


On Nov 8, 2009, at 3:11 PM, frenchcr wrote:




why do you use 365.25?



As opposed to what?

-- David


dates<-as.character(data[,"date_commissioned"]); # convert dates to
characters
#dates[1:10]
#[1] "19910101" "19860101" "19910101" "19860101" "19910101"  
"19910101"

"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
#dateObs[1:10]
#[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01"  
"1991-01-01"

"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"

today <- Sys.Date()
x.date <- as.Date(dateObs, format="%Y%m%d")

AGE <- round(as.vector(difftime(today , x.date, units='day') /  
365.25))






frenchcr wrote:



it sure does thank you!



will this work for you

x <- c('19910101', '19950302', '20010502')
today <- Sys.Date()
x.date <- as.Date(x, format="%Y%m%d")
round(as.vector(difftime(today , x.date, units='day') / 365.25))

[1] 19 15  9





On Sun, Nov 8, 2009 at 2:44 PM,   wrote:

Hi Jim,

Thanks for the quick reply...not sure what you mean by frame of
reference(only been using R for 4 days)...to clarify, i need to
turn my
dates from 1999-10-01 into 1999 then i subtract 2009 -1999 to get
an age
of 10. The column im working on has 312,000 rows and some have NA
in them
as we have no dates for that item.

To recap, the column is just a bunch of dates with some field  
empty, i
want to change the column from "date of commision" to "age of  
asset"


Cheers
Chris.





jholtman wrote:


What is the frame of reference to determine the age?   Check out
'difftime'.

On Sun, Nov 8, 2009 at 1:50 PM, frenchcr 
wrote:


Ive got a big column of dates (also some fields dont have a date
so they
have
NA instead),
that i have converted into date format as so...


dates<-as.character(data[,"date_commissioned"]); # converted  
dates to

characters
dates[1:10]
[1] "19910101" "19860101" "19910101" "19860101" "19910101"  
"19910101"

"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
dateObs[1:10]
[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01"  
"1991-01-01"

"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"



Now i need to turn the dates into AGE, how do i do it? Im not  
worried

about
fractions of years, whole years would do.


--
View this message in context:
http://old.nabble.com/Turn-dates-into-age- 
tp26256656p26256656.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.






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

What is the problem that 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.







--
View this message in context:
http://old.nabble.com/Turn-dates-into-age-tp26256656p26257435.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.


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

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

Re: [R] Simple 2-Way Anova issue in R

2009-11-08 Thread Jeremy Miles
If I've understood correctly, you have cell sizes of 1.  This is not enough.

ANOVA compares within group variance to between group variance, and
your within group variances are zero.

You need more data, or to collapse some cells.

Jeremy




2009/11/8 znd :
>
> Hello, I'm new to R and have been following many guides including the two-way
> anova (http://www.personality-project.org/r/r.anova.html).  Using that
> walkthrough including the supplied data I do str(data.ex2) and receive the
> appropriate types of data as follows:
>> str(data.ex2)
> 'data.frame':   16 obs. of  4 variables:
>  $ Observation: int  1 2 3 4 5 6 7 8 9 10 ...
>  $ Gender     : Factor w/ 2 levels "f","m": 2 2 2 2 2 2 2 2 1 1 ...
>  $ Dosage     : Factor w/ 2 levels "a","b": 1 1 1 1 2 2 2 2 1 1 ...
>  $ Alertness  : int  8 12 13 12 6 7 23 14 15 12 ...
>
> aov.ex2 = aov(Alertness~Gender*Dosage,data=data.ex2)
>
> summary(aov.ex2)
>
> Outputs:
>              Df  Sum Sq Mean Sq F value Pr(>F)
> Gender         1  76.562  76.562  2.9518 0.1115
> Dosage         1   5.062   5.062  0.1952 0.6665
> Gender:Dosage  1   0.063   0.063  0.0024 0.9617
> Residuals     12 311.250  25.938
>
> However, when I got to use my data that I made in csv format I have to tell
> R to interpret my factors which are year and depth as factors...
> datafilename="C:/Rclass/hmwk1pt2.csv"
> data.ex2=read.csv(datafilename,header=T)
> data.ex2$Year<-as.factor(data.ex2$Year)
> data.ex2$Depth<-as.factor(data.ex2$Depth)
> data.ex2
> str(data.ex2)
>
> This outputs what I would expect:
>
>> str(data.ex2)
> 'data.frame':   12 obs. of  4 variables:
>  $ Year      : Factor w/ 3 levels "1999","2000",..: 1 1 1 1 2 2 2 2 3 3 ...
>  $ Depth     : Factor w/ 4 levels "10","15","20",..: 1 2 3 4 1 2 3 4 1 2 ...
>  $ Replicate1: num  14.3 15.1 16.7 17.3 16.3 17.4 18.6 20.9 22.9 23.9 ...
>  $ Replicate2: num  14.7 15.6 16.9 17.9 16.4 17.2 19.6 21.3 22.7 23.3 ...
>
> But something is not causing my anova to carry through...this is what I
> have.
>
> ANOVA = aov(Replicate1~Year*Depth,data=data.ex2)
> summary(ANOVA)
>
> which outputs:
>
>> summary(ANOVA)
>            Df  Sum Sq Mean Sq
> Year         2 143.607  71.803
> Depth        3  17.323   5.774
> Year:Depth   6   2.587   0.431
>
> There is no F-value or Pr(>F) columns.
>
> I also can't boxplot this correctly, again following the example at that
> website above they have:
>
> boxplot(Alertness~Dosage*Gender,data=data.ex2)
>
> which outputs:
>
> http://old.nabble.com/file/p26258684/87o3uicpf6dt4kkdyvfv.jpeg
>
> My code is:
>
> boxplot(Replicate1~Year*Depth,data=data.ex2)
>
> which outputs:
>
> http://old.nabble.com/file/p26258684/gik02vyhvvbmcvw3ia2h.jpeg
>
> This is incorrect, it's multiplying my factors but I thought that when I did
> the str() on my data it recognized the Year and Depth as factors, not
> numbers or integers.
>
> My csv file is:
> http://old.nabble.com/file/p26258684/hmwk1pt2.csv hmwk1pt2.csv
>
> Any help on what is going one would be greatly appreciated because I need to
> perform one-way, two-way, nested, and factorial anovas but I first need to
> solve this problem before I can continue.
>
>
> --
> View this message in context: 
> http://old.nabble.com/Simple-2-Way-Anova-issue-in-R-tp26258684p26258684.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.
>



-- 
Jeremy Miles
Psychology Research Methods Wiki: www.researchmethodsinpsychology.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] ordered factor and unordered factor

2009-11-08 Thread Peng Yu
I don't understand under what situation ordered factor rather than
unordered factor should be used. Could somebody give me some examples?
What are the implications of order vs. unordered factors? Could
somebody recommend a textbook to me?

__
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] Turn dates into age

2009-11-08 Thread Jim Burke
To clarify.

Lets turn a date into an age. Given 05/29/1971 in mm/dd/
format. What is the year difference between then and today?
This would be the "age" requested that starts 05/29/1971 as
one.

Thanks,
Jim




David Winsemius wrote:
>
> On Nov 8, 2009, at 3:11 PM, frenchcr wrote:
>
>>
>>
>> why do you use 365.25?
>>
>
> As opposed to what?
>
> -- David
>>
>> dates<-as.character(data[,"date_commissioned"]); # convert dates to
>> characters
>> #dates[1:10]
>> #[1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
>> "19910101" "19910101" "19910101" "19910101"
>>
>> dateObs <- as.Date(dates,format="%Y%m%d")
>> #dateObs[1:10]
>> #[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
>> "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"
>>
>> today <- Sys.Date()
>> x.date <- as.Date(dateObs, format="%Y%m%d")
>>
>> AGE <- round(as.vector(difftime(today , x.date, units='day') / 365.25))
>>
>>
>>
>>
>>
>> frenchcr wrote:
>>>
>>>
>>> it sure does thank you!
>>>
>>>
 will this work for you

 x <- c('19910101', '19950302', '20010502')
 today <- Sys.Date()
 x.date <- as.Date(x, format="%Y%m%d")
 round(as.vector(difftime(today , x.date, units='day') / 365.25))
>>> [1] 19 15  9

>>>
>>>
>>> On Sun, Nov 8, 2009 at 2:44 PM,   wrote:
 Hi Jim,

 Thanks for the quick reply...not sure what you mean by frame of
 reference(only been using R for 4 days)...to clarify, i need to 
 turn my
 dates from 1999-10-01 into 1999 then i subtract 2009 -1999 to get 
 an age
 of 10. The column im working on has 312,000 rows and some have NA 
 in them
 as we have no dates for that item.

 To recap, the column is just a bunch of dates with some field empty, i
 want to change the column from "date of commision" to "age of asset"

 Cheers
 Chris.
>>>
>>>
>>>
>>>
>>> jholtman wrote:

 What is the frame of reference to determine the age?   Check out
 'difftime'.

 On Sun, Nov 8, 2009 at 1:50 PM, frenchcr  
 wrote:
>
> Ive got a big column of dates (also some fields dont have a date 
> so they
> have
> NA instead),
> that i have converted into date format as so...
>
>
> dates<-as.character(data[,"date_commissioned"]); # converted dates to
> characters
> dates[1:10]
> [1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
> "19910101" "19910101" "19910101" "19910101"
>
> dateObs <- as.Date(dates,format="%Y%m%d")
> dateObs[1:10]
> [1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
> "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"
>
>
>
> Now i need to turn the dates into AGE, how do i do it? Im not worried
> about
> fractions of years, whole years would do.
>
>
> -- 
> View this message in context:
> http://old.nabble.com/Turn-dates-into-age-tp26256656p26256656.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.
>



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

 What is the problem that 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.


>>>
>>>
>>
>> -- 
>> View this message in context: 
>> http://old.nabble.com/Turn-dates-into-age-tp26256656p26257435.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.
>
> David Winsemius, MD
> Heritage Laboratories
> West Hartford, CT
>
> __
> 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] Turn dates into age

2009-11-08 Thread David Winsemius


On Nov 8, 2009, at 3:11 PM, frenchcr wrote:




why do you use 365.25?



As opposed to what?

--  
David


dates<-as.character(data[,"date_commissioned"]); # convert dates to
characters
#dates[1:10]
#[1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
#dateObs[1:10]
#[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"

today <- Sys.Date()
x.date <- as.Date(dateObs, format="%Y%m%d")

AGE <- round(as.vector(difftime(today , x.date, units='day') /  
365.25))






frenchcr wrote:



it sure does thank you!



will this work for you

x <- c('19910101', '19950302', '20010502')
today <- Sys.Date()
x.date <- as.Date(x, format="%Y%m%d")
round(as.vector(difftime(today , x.date, units='day') / 365.25))

[1] 19 15  9





On Sun, Nov 8, 2009 at 2:44 PM,   wrote:

Hi Jim,

Thanks for the quick reply...not sure what you mean by frame of
reference(only been using R for 4 days)...to clarify, i need to  
turn my
dates from 1999-10-01 into 1999 then i subtract 2009 -1999 to get  
an age
of 10. The column im working on has 312,000 rows and some have NA  
in them

as we have no dates for that item.

To recap, the column is just a bunch of dates with some field  
empty, i

want to change the column from "date of commision" to "age of asset"

Cheers
Chris.





jholtman wrote:


What is the frame of reference to determine the age?   Check out
'difftime'.

On Sun, Nov 8, 2009 at 1:50 PM, frenchcr   
wrote:


Ive got a big column of dates (also some fields dont have a date  
so they

have
NA instead),
that i have converted into date format as so...


dates<-as.character(data[,"date_commissioned"]); # converted  
dates to

characters
dates[1:10]
[1] "19910101" "19860101" "19910101" "19860101" "19910101"  
"19910101"

"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
dateObs[1:10]
[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01"  
"1991-01-01"

"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"



Now i need to turn the dates into AGE, how do i do it? Im not  
worried

about
fractions of years, whole years would do.


--
View this message in context:
http://old.nabble.com/Turn-dates-into-age-tp26256656p26256656.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.





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

What is the problem that 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.







--
View this message in context: 
http://old.nabble.com/Turn-dates-into-age-tp26256656p26257435.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.


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] Simple 2-Way Anova issue in R

2009-11-08 Thread znd

Hello, I'm new to R and have been following many guides including the two-way
anova (http://www.personality-project.org/r/r.anova.html).  Using that
walkthrough including the supplied data I do str(data.ex2) and receive the
appropriate types of data as follows:
> str(data.ex2)
'data.frame':   16 obs. of  4 variables:
 $ Observation: int  1 2 3 4 5 6 7 8 9 10 ...
 $ Gender : Factor w/ 2 levels "f","m": 2 2 2 2 2 2 2 2 1 1 ...
 $ Dosage : Factor w/ 2 levels "a","b": 1 1 1 1 2 2 2 2 1 1 ...
 $ Alertness  : int  8 12 13 12 6 7 23 14 15 12 ...

aov.ex2 = aov(Alertness~Gender*Dosage,data=data.ex2)  

summary(aov.ex2)

Outputs:
  Df  Sum Sq Mean Sq F value Pr(>F)
Gender 1  76.562  76.562  2.9518 0.1115
Dosage 1   5.062   5.062  0.1952 0.6665
Gender:Dosage  1   0.063   0.063  0.0024 0.9617
Residuals 12 311.250  25.938 

However, when I got to use my data that I made in csv format I have to tell
R to interpret my factors which are year and depth as factors...
datafilename="C:/Rclass/hmwk1pt2.csv"
data.ex2=read.csv(datafilename,header=T)  
data.ex2$Year<-as.factor(data.ex2$Year)
data.ex2$Depth<-as.factor(data.ex2$Depth)
data.ex2
str(data.ex2)

This outputs what I would expect:

> str(data.ex2)
'data.frame':   12 obs. of  4 variables:
 $ Year  : Factor w/ 3 levels "1999","2000",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ Depth : Factor w/ 4 levels "10","15","20",..: 1 2 3 4 1 2 3 4 1 2 ...
 $ Replicate1: num  14.3 15.1 16.7 17.3 16.3 17.4 18.6 20.9 22.9 23.9 ...
 $ Replicate2: num  14.7 15.6 16.9 17.9 16.4 17.2 19.6 21.3 22.7 23.3 ...

But something is not causing my anova to carry through...this is what I
have.

ANOVA = aov(Replicate1~Year*Depth,data=data.ex2)
summary(ANOVA)

which outputs:

> summary(ANOVA)  
Df  Sum Sq Mean Sq
Year 2 143.607  71.803
Depth3  17.323   5.774
Year:Depth   6   2.587   0.431

There is no F-value or Pr(>F) columns.

I also can't boxplot this correctly, again following the example at that
website above they have:

boxplot(Alertness~Dosage*Gender,data=data.ex2)

which outputs:

http://old.nabble.com/file/p26258684/87o3uicpf6dt4kkdyvfv.jpeg 

My code is:

boxplot(Replicate1~Year*Depth,data=data.ex2)

which outputs:

http://old.nabble.com/file/p26258684/gik02vyhvvbmcvw3ia2h.jpeg 

This is incorrect, it's multiplying my factors but I thought that when I did
the str() on my data it recognized the Year and Depth as factors, not
numbers or integers.

My csv file is:
http://old.nabble.com/file/p26258684/hmwk1pt2.csv hmwk1pt2.csv 

Any help on what is going one would be greatly appreciated because I need to
perform one-way, two-way, nested, and factorial anovas but I first need to
solve this problem before I can continue.


-- 
View this message in context: 
http://old.nabble.com/Simple-2-Way-Anova-issue-in-R-tp26258684p26258684.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] rd doc truncated with R 2.10.0

2009-11-08 Thread Patrick Giraudoux

Patrick Giraudoux a écrit :

Duncan Murdoch a écrit :

On 08/11/2009 12:07 PM, Patrick Giraudoux wrote:

Hi,

I am routinely compiling a package and since I have moved to R 
2.10.0, it troncates some section texts in the doc:


With the following section in the rd file:

\details{
 The function calls gpsbabel via the system. The gpsbabel program 
must be present and on the user's PATH for the function to work see 
. The function has been tested on the 
following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and 
GPSmap 60CSx.

}

...compiling under R 2.9.2 (rcmd build --binary --auto-zip pgirmess) 
I get this


Details:

 The function calls gpsbabel via the system. The gpsbabel program
 must be present and on the user's PATH for the function to work,
 see . The function has been tested on
 the following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and
 GPSmap 60CSx.

and compiling now under R 2.10.0

Details:

 The function has been tested on the following Garmin GPS devices:
 Etrex Summit, Etrex Vista Cx and GPSmap 60CSx. The function calls
 gpsbabel via the system. The gpsbabel program must be presen


Has anyone an explanation and a workaround ?


You will need to make the complete file available to us to diagnose 
this.  Is it in pgirmess 1.4.0?  Which topic?


Duncan Murdoch


OK. Will send it offlist.

Patrick



Just to end up the thread, here is the solution, thanks to Duncan Murdoch:

"It was an embarassing bug in the C code that did tab expansion:  it 
could only handle lines with 200 chars in them, and yours had about 270, 
so the last part was lost.


It will soon be fixed in R-patched, to become 2.10.1 in due course.

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.


Re: [R] reference on contr.helmert and typo on its help page.

2009-11-08 Thread Peng Yu
On Sun, Nov 8, 2009 at 11:28 AM, Peter Dalgaard
 wrote:
> Gabor Grothendieck wrote:
>>
>> On Sun, Nov 8, 2009 at 11:59 AM, Peng Yu  wrote:
>>>
>>> On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch 
>>> wrote:

 On 08/11/2009 11:03 AM, Peng Yu wrote:
>
> I'm wondering which textbook discussed the various contrast matrices
> mentioned in the help page of 'contr.helmert'. Could somebody let me
> know?

 Doesn't the reference on that page discuss them?
>>>
>>> It does explain what the functions are. But I need a more basic and
>>> complete reference. For example, I want to understand what 'Helmert
>>> parametrization' (on page 33 of 'Statistical Models in S') is.
>>>
>>
>> Just google for: Helmert contrasts
>
> Or,
>
>> contr.helmert(5)
>  [,1] [,2] [,3] [,4]
> 1   -1   -1   -1   -1
> 2    1   -1   -1   -1
> 3    0    2   -1   -1
> 4    0    0    3   -1
> 5    0    0    0    4
>
>> MASS::fractions(MASS::ginv(contr.helmert(5)))
>     [,1]  [,2]  [,3]  [,4]  [,5]
> [1,]  -1/2   1/2     0     0     0
> [2,]  -1/6  -1/6   1/3     0     0
> [3,] -1/12 -1/12 -1/12   1/4     0
> [4,] -1/20 -1/20 -1/20 -1/20   1/5
>
> and apply brains.
>
> I.e., except for a slightly odd multiplier, the parameters represent the
>  difference between each level and the average of the preceding levels.

I realized that my questions are what a contrast matrix is and how it
is related to hypothesis testing. For a give hypothesis, how to get
the corresponding contrast matrix in a systematical way? There are
some online materials, but they are all diffused. I have also read the
book Applied Linear Regression Models, which doesn't give a complete
descriptions on all the aspects of contrast and contrast matrix. But I
would want a textbook that gives a complete description, so that I
don't have to look around for other materials.

__
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] Obtaining midpoints of class intervals produced by cut and table

2009-11-08 Thread Gabor Grothendieck
Try this:

library(gsubfn)
demo("gsubfn-cut")

and note the strapply call that converts the levels from cut to a matrix.

On Sun, Nov 8, 2009 at 4:08 PM, jose romero  wrote:
> Hello list:
>
> I am using "cut" and "table" to obtain a frequency table from a numeric 
> sample vector.  The idea is to calculate mean and standard deviation on 
> grouped data.  However, I can't extract the midpoints of the class intervals, 
> which seem to be strings treated as factors.  How do i extract the midpoint?
>
> Thanks,
> jose loreto
>
>
>        [[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] Obtaining midpoints of class intervals produced by cut and table

2009-11-08 Thread Peter Dalgaard

jose romero wrote:

Hello list:

I am using "cut" and "table" to obtain a frequency table from a

numeric sample vector. The idea is to calculate mean and standard
deviation on grouped data. However, I can't extract the midpoints of the
class intervals, which seem to be strings treated as factors. How do i
extract the midpoint?




You might consider starting with hist() instead. Otherwise, how about

mid <- brk[-1] - diff(brk)/2

where brk is the breakpoints used in cut()

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

__
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] Obtaining midpoints of class intervals produced by cut and table

2009-11-08 Thread baptiste auguie
Hi,

Maybe something like this (inspired by ?cut),

cut2num <- function(f){
  labs <- levels(f)
  d <- data.frame(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
  upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) ))
  d$midpoints <- rowMeans(d)
  d
}

a <- c(1, 2, 3, 4, 5, 2, 3, 4, 5, 6, 7)
cut2num(cut(a, 3))

HTH,

baptiste





2009/11/8 jose romero :
> Hello list:
>
> I am using "cut" and "table" to obtain a frequency table from a numeric 
> sample vector.  The idea is to calculate mean and standard deviation on 
> grouped data.  However, I can't extract the midpoints of the class intervals, 
> which seem to be strings treated as factors.  How do i extract the midpoint?
>
> Thanks,
> jose loreto
>
>
>        [[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] Obtaining midpoints of class intervals produced by cut and table

2009-11-08 Thread jose romero
Hello list:

I am using "cut" and "table" to obtain a frequency table from a numeric sample 
vector.  The idea is to calculate mean and standard deviation on grouped data.  
However, I can't extract the midpoints of the class intervals, which seem to be 
strings treated as factors.  How do i extract the midpoint?

Thanks,
jose loreto


[[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] Turn dates into age

2009-11-08 Thread frenchcr


why do you use 365.25?


dates<-as.character(data[,"date_commissioned"]); # convert dates to
characters
#dates[1:10]
#[1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
#dateObs[1:10]
#[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"

today <- Sys.Date()
x.date <- as.Date(dateObs, format="%Y%m%d")

AGE <- round(as.vector(difftime(today , x.date, units='day') / 365.25))





frenchcr wrote:
> 
> 
> it sure does thank you!
> 
> 
>> will this work for you
>>
>> x <- c('19910101', '19950302', '20010502')
>> today <- Sys.Date()
>> x.date <- as.Date(x, format="%Y%m%d")
>> round(as.vector(difftime(today , x.date, units='day') / 365.25))
> [1] 19 15  9
>>
> 
> 
> On Sun, Nov 8, 2009 at 2:44 PM,   wrote:
>> Hi Jim,
>>
>> Thanks for the quick reply...not sure what you mean by frame of
>> reference(only been using R for 4 days)...to clarify, i need to turn my
>> dates from 1999-10-01 into 1999 then i subtract 2009 -1999 to get an age
>> of 10. The column im working on has 312,000 rows and some have NA in them
>> as we have no dates for that item.
>>
>> To recap, the column is just a bunch of dates with some field empty, i
>> want to change the column from "date of commision" to "age of asset"
>>
>> Cheers
>> Chris.
> 
> 
> 
> 
> jholtman wrote:
>> 
>> What is the frame of reference to determine the age?   Check out
>> 'difftime'.
>> 
>> On Sun, Nov 8, 2009 at 1:50 PM, frenchcr  wrote:
>>>
>>> Ive got a big column of dates (also some fields dont have a date so they
>>> have
>>> NA instead),
>>> that i have converted into date format as so...
>>>
>>>
>>> dates<-as.character(data[,"date_commissioned"]); # converted dates to
>>> characters
>>> dates[1:10]
>>> [1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
>>> "19910101" "19910101" "19910101" "19910101"
>>>
>>> dateObs <- as.Date(dates,format="%Y%m%d")
>>> dateObs[1:10]
>>> [1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
>>> "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"
>>>
>>>
>>>
>>> Now i need to turn the dates into AGE, how do i do it? Im not worried
>>> about
>>> fractions of years, whole years would do.
>>>
>>>
>>> --
>>> View this message in context:
>>> http://old.nabble.com/Turn-dates-into-age-tp26256656p26256656.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.
>>>
>> 
>> 
>> 
>> -- 
>> Jim Holtman
>> Cincinnati, OH
>> +1 513 646 9390
>> 
>> What is the problem that 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.
>> 
>> 
> 
> 

-- 
View this message in context: 
http://old.nabble.com/Turn-dates-into-age-tp26256656p26257435.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] Turn dates into age

2009-11-08 Thread frenchcr


it sure does thank you!


> will this work for you
>
> x <- c('19910101', '19950302', '20010502')
> today <- Sys.Date()
> x.date <- as.Date(x, format="%Y%m%d")
> round(as.vector(difftime(today , x.date, units='day') / 365.25))
[1] 19 15  9
>


On Sun, Nov 8, 2009 at 2:44 PM,   wrote:
> Hi Jim,
>
> Thanks for the quick reply...not sure what you mean by frame of
> reference(only been using R for 4 days)...to clarify, i need to turn my
> dates from 1999-10-01 into 1999 then i subtract 2009 -1999 to get an age
> of 10. The column im working on has 312,000 rows and some have NA in them
> as we have no dates for that item.
>
> To recap, the column is just a bunch of dates with some field empty, i
> want to change the column from "date of commision" to "age of asset"
>
> Cheers
> Chris.




jholtman wrote:
> 
> What is the frame of reference to determine the age?   Check out
> 'difftime'.
> 
> On Sun, Nov 8, 2009 at 1:50 PM, frenchcr  wrote:
>>
>> Ive got a big column of dates (also some fields dont have a date so they
>> have
>> NA instead),
>> that i have converted into date format as so...
>>
>>
>> dates<-as.character(data[,"date_commissioned"]); # converted dates to
>> characters
>> dates[1:10]
>> [1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
>> "19910101" "19910101" "19910101" "19910101"
>>
>> dateObs <- as.Date(dates,format="%Y%m%d")
>> dateObs[1:10]
>> [1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
>> "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"
>>
>>
>>
>> Now i need to turn the dates into AGE, how do i do it? Im not worried
>> about
>> fractions of years, whole years would do.
>>
>>
>> --
>> View this message in context:
>> http://old.nabble.com/Turn-dates-into-age-tp26256656p26256656.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.
>>
> 
> 
> 
> -- 
> Jim Holtman
> Cincinnati, OH
> +1 513 646 9390
> 
> What is the problem that 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.
> 
> 

-- 
View this message in context: 
http://old.nabble.com/Turn-dates-into-age-tp26256656p26257419.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] MCMC gradually slows down

2009-11-08 Thread Jens Malmros
Thanks a lot, this works.

jim holtman wrote:
> First of all, allocate 'theta' to be the final size you need.  Every
> time through your loop you are extending it by one, meaning you are
> spending a lot of time copying the data each time.  Do something like:
> 
> theta <- numeric(n)
> 
> and then see how fast it works.
> 
> On Sun, Nov 8, 2009 at 2:11 PM, Jens Malmros  wrote:
>> Hello,
>>
>> I have written a simple Metropolis-Hastings MCMC algorithm for a
>> binomial parameter:
>>
>> MHastings = function(n,p0,d){
>>theta = c()
>>theta[1] = p0
>>t =1
>>while(t<=n){
>>phi = log(theta[t]/(1-theta[t]))
>>phisim = phi + rnorm(1,0,d)
>>thetasim = exp(phisim)/(1+exp(phisim))
>>r = (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8)
>>if(runif(1,0,1)>theta[t+1] = thetasim
>>} else {
>>theta[t+1] = theta[t]
>>}
>>t = t+1
>>if(t%%1000==0) print(t) # diagnostic
>>}
>>data.frame(theta)
>> }
>>
>> The problem is that it gradually slows down. It is very fast in the
>> beginning, but slows down and gets very slow as you reach about 5
>> iterations and I need do to plenty more.
>>
>> I know there are more fancy MCMC routines available, but I am really
>> just interested in this to work.
>>
>> Thank you for your help,
>> Jens Malmros
>>
>> __
>> 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] Models for Discrete Choice in R

2009-11-08 Thread Emmanuel Charpentier
Le dimanche 08 novembre 2009 à 17:07 -0200, Iuri Gavronski a écrit :
> Hi,
> 
> I would like to fit Logit models for ordered data, such as those
> suggested by Greene (2003), p. 736.
> 
> Does anyone suggests any package in R for that?

look up the polr function in package MASS (and read the relevant pages
in V&R4 and some quoted references...) or the slightly more
sophisticated (larger range of models) lrm function in F. Harrell's
Design (now rms) packge (but be aware that Design is a huge beast witch
carries its own "computing universe", based on (strong) Harrell's view
of what a regression analysis should be : reading his book is, IMHO,
necessary to understand his choices and agree (or disgree) with them).

If you have a multilevel model (a. k. a. one "random effect" grouping),
the "repolr" packge aims at that, but I've been unable to use it
recently (numerical exceptions).

> By the way, my dependent variable is ordinal and my independent
> variables are ratio/intervalar.

Numeric ? Then maybe some recoding/transformation is in order ... in
which case Design/rms might or might not be useful.

HTH,

Emmanuel Charpentier

__
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] MCMC gradually slows down

2009-11-08 Thread William Dunlap
> -Original Message-
> From: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] On Behalf Of Jens Malmros
> Sent: Sunday, November 08, 2009 11:11 AM
> To: r-help@r-project.org
> Subject: [R] MCMC gradually slows down
> 
> Hello,
> 
> I have written a simple Metropolis-Hastings MCMC algorithm for a
> binomial parameter:
> 
> MHastings = function(n,p0,d){
>   theta = c()

Change that theta <- c() to
  theta <- numeric(n)
and it will go faster.  Growing datasets
can result in a lot of unnecessary memory
management.  The original had an obvious
quadratic quality in the plot of n vs. MHastings(n,.2,2)
and preallocating theta to its final length
made it look linear up to n=10 and at n=10
the time for the original was 35 seconds vs 3.75
for the preallocated version.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

>   theta[1] = p0
>   t =1
>   while(t<=n){
>   phi = log(theta[t]/(1-theta[t]))
>   phisim = phi + rnorm(1,0,d)
>   thetasim = exp(phisim)/(1+exp(phisim))
>   r = 
> (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8)
>   if(runif(1,0,1)   theta[t+1] = thetasim
>   } else {
>   theta[t+1] = theta[t]
>   }
>   t = t+1
>   if(t%%1000==0) print(t) # diagnostic
>   }
>   data.frame(theta)
> }
> 
> The problem is that it gradually slows down. It is very fast in the
> beginning, but slows down and gets very slow as you reach about 5
> iterations and I need do to plenty more.
> 
> I know there are more fancy MCMC routines available, but I am really
> just interested in this to work.
> 
> Thank you for your help,
> Jens Malmros
> 
> __
> 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] a question in coxph

2009-11-08 Thread Lei Liu

Hi there,

I tried to fit a penalized spline for a continuous risk factor in 
recurrent events data. I found that I can include both pspline and 
frailty terms in coxph. So I use code like


  fit1 <- coxph(Surv(start, end, event)~pspline(age, df=0) + male + 
white +frailty(id), data.all)


It yields results for both age spline and frailty variance. Next I 
tried to draw the plot for age effect. I used termplot function


  termplot(fit1, terms=c(1), se=TRUE)

However, it gives me some error message:

Error in dimnames(pred) <- list(dimnames(x)[[1]], termname) :
  length of 'dimnames' [2] not equal to array extent

I don't know why. However, if I delete the frailty option, and use only

  fit2 <- coxph(Surv(start, end, event)~pspline(age, df=0) + male + 
white, data.all)


I can generate the plot by

termplot(fit2, terms=c(1), se=TRUE)

So it seems to me that adding the frailty() option makes the 
"termplot" function fail to work? Do you know how to draw the plot 
for my function fit1? Thanks a lot!


Lei


At 09:53 AM 7/20/2009, you wrote:

 I would be willing to chair the session.
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.


Re: [R] Turn dates into age

2009-11-08 Thread Marc Schwartz
As Jim has noted, if the dates you have below are an 'end date', you  
need to define the time0 or start date for each to calculate the  
intervals. On the other hand, are the dates you have below the start  
dates and you need to calculate the time to today? In the latter case,  
see ?Sys.Date to get the current system date from your computer.


Generically speaking you can use the following:

(EndDate - StartDate) / 365.25

where EndDate and StartDate are of class Date. This will give you the  
time interval in years plus any fraction.


You can then use round() which will give you a whole year with typical  
rounding up or down to the nearest whole integer. You can use floor(),  
which will give you the nearest whole integer less than the result or  
basically a round down result. Keep in mind that the above calculation  
does not honor a calendar year, but is an approximation.


If you want to calculate age in years as we typically think of it  
using the calendar, you can use the following, where DOB is the Date  
of Birth (Start Date) and Date2 is the End Date:


# Calculate Age in Years
# DOB: Class Date
# Date2: Class Date

Calc_Age <- function(DOB, Date2)
{
  if (length(DOB) != length(Date2))
stop("length(DOB) != length(Date2)")

  if (!inherits(DOB, "Date") | !inherits(Date2, "Date"))
stop("Both DOB and Date2 must be Date class objects")

  start <- as.POSIXlt(DOB)
  end <- as.POSIXlt(Date2)

  Age <- end$year - start$year

  ifelse((end$mon < start$mon) |
 ((end$mon == start$mon) & (end$mday < start$mday)),
 Age - 1, Age)
}


HTH,

Marc Schwartz


On Nov 8, 2009, at 1:30 PM, jim holtman wrote:

What is the frame of reference to determine the age?   Check out  
'difftime'.


On Sun, Nov 8, 2009 at 1:50 PM, frenchcr   
wrote:


Ive got a big column of dates (also some fields dont have a date so  
they have

NA instead),
that i have converted into date format as so...


dates<-as.character(data[,"date_commissioned"]); # converted dates to
characters
dates[1:10]
[1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
dateObs[1:10]
[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"



Now i need to turn the dates into AGE, how do i do it? Im not  
worried about

fractions of years, whole years would do.


__
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] Extracting matched expressions

2009-11-08 Thread Gabor Grothendieck
strapply in the gsubfn package can do that. It applies the indicated
function, here just c, to the back references from the pattern match
and then simplifies the result using simplify. (If you omit simplify
here it would give a one element list like strsplit does.)

library(gsubfn)
pat <- "(.*?) (.*?) ([ehtr]{5})"
strapply("one two three", pat, c, simplify = c)

See home page at: http://gsubfn.googlecode.com


On Sun, Nov 8, 2009 at 1:51 PM, Hadley Wickham  wrote:
> Hi all,
>
> Is there a tool in base R to extract matched expressions from a
> regular expression?  i.e. given the regular expression "(.*?) (.*?)
> ([ehtr]{5})" is there a way to extract the character vector c("one",
> "two", "three") from the string "one two three" ?
>
> Thanks,
>
> 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.
>

__
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] MCMC gradually slows down

2009-11-08 Thread Viechtbauer Wolfgang (STAT)
Try this:

Instead of:

theta = c()

use:

theta <- rep(NA, 50)

or however many iterations you want the algorithm to run.

Best,

--
Wolfgang Viechtbauerhttp://www.wvbauer.com/
Department of Methodology and StatisticsTel: +31 (0)43 388-2277
School for Public Health and Primary Care   Office Location:
Maastricht University, P.O. Box 616 Room B2.01 (second floor)
6200 MD Maastricht, The Netherlands Debyeplein 1 (Randwyck)

From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf Of 
Jens Malmros [jens.malm...@gmail.com]
Sent: Sunday, November 08, 2009 8:11 PM
To: r-help@r-project.org
Subject: [R] MCMC gradually slows down

Hello,

I have written a simple Metropolis-Hastings MCMC algorithm for a
binomial parameter:

MHastings = function(n,p0,d){
theta = c()
theta[1] = p0
t =1
while(t<=n){
phi = log(theta[t]/(1-theta[t]))
phisim = phi + rnorm(1,0,d)
thetasim = exp(phisim)/(1+exp(phisim))
r = (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8)
if(runif(1,0,1)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] MCMC gradually slows down

2009-11-08 Thread jim holtman
First of all, allocate 'theta' to be the final size you need.  Every
time through your loop you are extending it by one, meaning you are
spending a lot of time copying the data each time.  Do something like:

theta <- numeric(n)

and then see how fast it works.

On Sun, Nov 8, 2009 at 2:11 PM, Jens Malmros  wrote:
> Hello,
>
> I have written a simple Metropolis-Hastings MCMC algorithm for a
> binomial parameter:
>
> MHastings = function(n,p0,d){
>        theta = c()
>        theta[1] = p0
>        t =1
>        while(t<=n){
>                phi = log(theta[t]/(1-theta[t]))
>                phisim = phi + rnorm(1,0,d)
>                thetasim = exp(phisim)/(1+exp(phisim))
>                r = (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8)
>                if(runif(1,0,1)                        theta[t+1] = thetasim
>                } else {
>                        theta[t+1] = theta[t]
>                }
>                t = t+1
>                if(t%%1000==0) print(t) # diagnostic
>        }
>        data.frame(theta)
> }
>
> The problem is that it gradually slows down. It is very fast in the
> beginning, but slows down and gets very slow as you reach about 5
> iterations and I need do to plenty more.
>
> I know there are more fancy MCMC routines available, but I am really
> just interested in this to work.
>
> Thank you for your help,
> Jens Malmros
>
> __
> 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 that 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.


Re: [R] Extracting matched expressions

2009-11-08 Thread jim holtman
Is this what you want:

> x <- '  one two three  '
> y <- 
> sub(".*?([^[:space:]]+)[[:space:]]+([^[:space:]]+)[[:space:]]+([ehrt]{5}).*",
+ "\\1 \\2 \\3", x, perl=TRUE)
> unlist(strsplit(y, ' '))
[1] "one"   "two"   "three"


On Sun, Nov 8, 2009 at 1:51 PM, Hadley Wickham  wrote:
> Hi all,
>
> Is there a tool in base R to extract matched expressions from a
> regular expression?  i.e. given the regular expression "(.*?) (.*?)
> ([ehtr]{5})" is there a way to extract the character vector c("one",
> "two", "three") from the string "one two three" ?
>
> Thanks,
>
> 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.
>



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

What is the problem that 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.


Re: [R] correlated binary data and overall probability

2009-11-08 Thread Christian Lerch

To answer my own question:
The package mvtBinaryEP was helpful!

Christian Lerch schrieb:

Dear All,

I try to simulate correlated binary data for a clinical research project.
Unfortunately, I do not come to grips with bindata().

Consider
corr<-data.frame(ID=as.factor(rep(c(1:10), each=5)),
task=as.factor(rep(c(1:5),10)))

[this format might be more appropriate:]
corr2<-data.frame(ID=as.factor(rep(c(1:10),5)),
tablet=as.factor(rep(c(1:5),each=10)))

Now, I want to add one column 'outcome' for the binary response:
* within each 'task', the probabilty of success (1) is fixed, (say 
rbinom(x,1,0.7))

* within each 'ID', the outcomes are correlated (say 0.8)

How can I generate this column 'outcome' with the given proporties?

Many thanks for hints or even code!

Regards,
Christian


__
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] Turn dates into age

2009-11-08 Thread jim holtman
What is the frame of reference to determine the age?   Check out 'difftime'.

On Sun, Nov 8, 2009 at 1:50 PM, frenchcr  wrote:
>
> Ive got a big column of dates (also some fields dont have a date so they have
> NA instead),
> that i have converted into date format as so...
>
>
> dates<-as.character(data[,"date_commissioned"]); # converted dates to
> characters
> dates[1:10]
> [1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
> "19910101" "19910101" "19910101" "19910101"
>
> dateObs <- as.Date(dates,format="%Y%m%d")
> dateObs[1:10]
> [1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
> "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"
>
>
>
> Now i need to turn the dates into AGE, how do i do it? Im not worried about
> fractions of years, whole years would do.
>
>
> --
> View this message in context: 
> http://old.nabble.com/Turn-dates-into-age-tp26256656p26256656.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.
>



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

What is the problem that 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.


Re: [R] lme4 and incomplete block design

2009-11-08 Thread Christian Lerch

Many thanks, Bill and Emmanuel!
Christian

Emmanuel Charpentier schrieb:

Le dimanche 08 novembre 2009 à 00:05 +0100, Christian Lerch a écrit :

Dear list members,

I try to simulate an incomplete block design in which every participants 
receives 3 out of 4 possible treatment. The outcome in binary.


Assigning a binary outcome to the BIB or PBIB dataset of the package 
SASmixed gives the appropriate output.
With the code below, fixed treatment estimates are not given for each of 
the 4 possible treatments, instead a kind of summary measure(?) for 
'treatment' is given.


block<-rep(1:24,each=3)
treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 
2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 
1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 
4, 4, 1, 3)
outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
0, 0, 1, 0)

data<-data.frame(block,treatment,outcome)
lmer(outcome~treatment +(1|block), family=binomial, data=data)

Is this a problem with the recovery of interblock estimates?


No...

  Are there 
special rules how incomplete block designs should look like to enable 
this recovery?


Neither...

Compare :


library(lme4)

Le chargement a nécessité le package : Matrix
Le chargement a nécessité le package : lattice

summary(lmer(outcome~treatment +(1|block), family=binomial,

+  data=data.frame(block<-rep(1:24,each=3),
+treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4,
+ 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4,
+ 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4,
+ 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2,
+ 4, 4, 1, 3),
+outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 
0,
+   0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 
0,
+   0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 
1,
+   0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 
0))
+  ))
Generalized linear mixed model fit by the Laplace approximation 
Formula: outcome ~ treatment + (1 | block) 
   Data: data.frame(block <- rep(1:24, each = 3), treatment <- c(1, 2,  3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2,  4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2,  1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2,  1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3), outcome <- c(0, 0,  0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,  0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0,  0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,  0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) 
   AIC  BIC logLik deviance

 86.28 93.1 -40.1480.28
Random effects:
 Groups NameVariance Std.Dev.
 block  (Intercept) 0.60425  0.77734 
Number of obs: 72, groups: block, 24


Fixed effects:
Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.277830.71767  -1.7800.075 .
treatment0.011620.25571   0.0450.964  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


Correlation of Fixed Effects:
  (Intr)
treatment -0.892

with :


summary(lmer(outcome~treatment +(1|block), family=binomial,

+  data=data.frame(block<-rep(1:24,each=3),
+treatment<-factor(c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 
4,
+2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 
1,
+3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 
1,
+4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 
1,
+1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)),
+outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 
0,
+   0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 
0,
+   0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 
1,
+   0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 
0))
+  ))
Generalized linear mixed model fit by the Laplace approximation 
Formula: outcome ~ treatment + (1 | block) 
   Data: data.frame(block <- rep(1:24, each = 3), treatment <- factor(c(1,  2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1,  2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3,  2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4,  2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)), outcome <- c(0,  0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,  0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1,  0, 0, 0, 0, 0, 0, 0, 1, 0

[R] Turn dates into age

2009-11-08 Thread frenchcr


Ive got a big column of dates (also some fields dont have a date so they
have NA instead),
that i have converted into date format as so...


dates<-as.character(data[,"date_commissioned"]); # converted dates to
characters
dates[1:10]
[1] "19910101" "19860101" "19910101" "19860101" "19910101" "19910101"
"19910101" "19910101" "19910101" "19910101"

dateObs <- as.Date(dates,format="%Y%m%d")
dateObs[1:10]
[1] "1991-01-01" "1986-01-01" "1991-01-01" "1986-01-01" "1991-01-01"
"1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01" "1991-01-01"



Now i need to turn the dates into years, how do i do it? Im not worried
about fractions of years, whole years would do.


-- 
View this message in context: 
http://old.nabble.com/Turn-dates-into-age-tp26256656p26256656.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] negative log likelihood

2009-11-08 Thread mat7770

I have two related variables, each with 16 points (x and Y). I am given
variance and the y-intercept. I know how to create a regression line and
find the residuals, but here is my problem. I have to make a loop that uses
the seq() function, so that it changes the slope value of the y=mx + B
equation ranging from 0-5 in increments of 0.01. The loop also needs to
calculate the negative log likelihood at each slope value and determine the
lowest one. I know that R can compute the best regression line by using
lm(y~x), but I need to see if that value matches the loop functions. 
-- 
View this message in context: 
http://old.nabble.com/negative-log-likelihood-tp26256881p26256881.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] MCMC gradually slows down

2009-11-08 Thread Jens Malmros
Hello,

I have written a simple Metropolis-Hastings MCMC algorithm for a
binomial parameter:

MHastings = function(n,p0,d){
theta = c()
theta[1] = p0
t =1
while(t<=n){
phi = log(theta[t]/(1-theta[t]))
phisim = phi + rnorm(1,0,d)
thetasim = exp(phisim)/(1+exp(phisim))
r = (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8)
if(runif(1,0,1)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] Models for Discrete Choice in R

2009-11-08 Thread Iuri Gavronski
Hi,

I would like to fit Logit models for ordered data, such as those
suggested by Greene (2003), p. 736.

Does anyone suggests any package in R for that?

By the way, my dependent variable is ordinal and my independent
variables are ratio/intervalar.

Thanks,

Iuri.

Greene, W. H. Econometric Analysis. Upper Saddle River, NJ: Prentice Hall, 2003

__
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] Extracting matched expressions

2009-11-08 Thread Hadley Wickham
Hi all,

Is there a tool in base R to extract matched expressions from a
regular expression?  i.e. given the regular expression "(.*?) (.*?)
([ehtr]{5})" is there a way to extract the character vector c("one",
"two", "three") from the string "one two three" ?

Thanks,

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] How to create multiple command windows of R on Windows?

2009-11-08 Thread Ivan
Thanks a lot for your explanation. That actually makes sense.

On Sat, Nov 7, 2009 at 3:54 PM, Duncan Murdoch  wrote:

> On 07/11/2009 6:31 PM, Ivan wrote:
>
>> Well, the problem is that I want those console to be from the same session
>> (i.e., they share same objects etc.). I do not think multiple instances of
>> R
>> could do it ...
>>
>
> No, there's only one console in each instance.  That's pretty fundamental:
>  there is one master loop reading, evaluating and printing, and the console
> shows you what's going on there.
>
> You can write your own front end to send things to the console (and somehow
> decide what to do when it's busy), but don't expect R to ever have more than
> one console.
>
> Duncan Murdoch
>
>
>  2009/11/7 Uwe Ligges 
>>
>>  You can start as many instances of R as you like (except for memory
>>> restrictions, perhaps), hence I do not get your point.
>>>
>>> Uwe Ligges
>>>
>>>
>>>
>>>
>>> Ivan wrote:
>>>
>>>  Thanks for that. I seem to be able to only get one R console here
 though.
 Actually that is my question: how to get a different R console?


 On Fri, Nov 6, 2009 at 4:10 PM,  wrote:

  It's easy.

> You can execute R commands under different ``command consoles'' in
> Windows.
>
>
> regards
>
>
>  Huang, Guo-Hao
>
> --
> From: "Ivan" 
> Sent: Saturday, November 07, 2009 8:01 AM
> To: 
> Subject: [R] How to create multiple command windows of R on Windows?
>
> Hi there,
>
>
>
> I am using R on Windows here. And I got a rather stupid question: how
> can
> I
> create another R console? It would be nice to have multiple R consoles
> so
> that I can separate different types of commands I use.
>
>
>
> Thank you,
>
>
>
> Ivan
>
> [[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.


>>[[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] rd doc truncated with R 2.10.0

2009-11-08 Thread Patrick Giraudoux

Duncan Murdoch a écrit :

On 08/11/2009 12:07 PM, Patrick Giraudoux wrote:

Hi,

I am routinely compiling a package and since I have moved to R 
2.10.0, it troncates some section texts in the doc:


With the following section in the rd file:

\details{
 The function calls gpsbabel via the system. The gpsbabel program 
must be present and on the user's PATH for the function to work see 
. The function has been tested on the 
following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and GPSmap 
60CSx.

}

...compiling under R 2.9.2 (rcmd build --binary --auto-zip pgirmess) 
I get this


Details:

 The function calls gpsbabel via the system. The gpsbabel program
 must be present and on the user's PATH for the function to work,
 see . The function has been tested on
 the following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and
 GPSmap 60CSx.

and compiling now under R 2.10.0

Details:

 The function has been tested on the following Garmin GPS devices:
 Etrex Summit, Etrex Vista Cx and GPSmap 60CSx. The function calls
 gpsbabel via the system. The gpsbabel program must be presen


Has anyone an explanation and a workaround ?


You will need to make the complete file available to us to diagnose 
this.  Is it in pgirmess 1.4.0?  Which topic?


Duncan Murdoch


OK. Will send it offlist.

Patrick

__
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] linear trend line and a quadratic trend line.

2009-11-08 Thread Peter Dalgaard

Eric Fail wrote:



Dear list users



How is it possible to visualise both a linear trend line and a quadratic trend 
line on a plot
of two variables?



Here my almost working exsample.



data(Duncan)

attach(Duncan)

plot(prestige ~ income)

abline(lm(prestige ~ income), col=2, lwd=2)



Now I would like to add yet another trend line, but this time a quadratic one. 
So I have two
trend lines. One linear trend line and a quadratic trend line. A trend line as 
if I had taken
I(income^2) into the equation.



I know I can make two models and compare them using anova, but for pedagogical 
resons I wold
like to visualise it. I know the trick from my past as an SPSS user, but I 
would like to do
this in R as well. Se it in SPSS
http://www.childrens-mercy.org/stats/weblog2006/QuadraticRegression.asp


There's no precooked function that I am aware of, but the generic way is 
like


rg <- range(income)
N <- 200
x <- seq(rg[1], rg[2],, N)
pred <- predict(lm(prestige~ income+I(income^2)),
  newdata=data.frame(income=x))
lines(x, pred)

as usual, "like" means that if you can't be bothered with making your 
example reproducible, I can't be bothered with testing the code!



Well, actually, I found the Duncan data in library(car), so I did in 
fact test...


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

__
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] rd doc truncated with R 2.10.0

2009-11-08 Thread Duncan Murdoch

On 08/11/2009 12:07 PM, Patrick Giraudoux wrote:

Hi,

I am routinely compiling a package and since I have moved to R 2.10.0, 
it troncates some section texts in the doc:


With the following section in the rd file:

\details{
 The function calls gpsbabel via the system. The gpsbabel program must 
be present and on the user's PATH for the function to work see 
. The function has been tested on the 
following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and GPSmap 60CSx.

}

...compiling under R 2.9.2 (rcmd build --binary --auto-zip pgirmess) I 
get this


Details:

 The function calls gpsbabel via the system. The gpsbabel program
 must be present and on the user's PATH for the function to work,
 see . The function has been tested on
 the following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and
 GPSmap 60CSx.

and compiling now under R 2.10.0

Details:

 The function has been tested on the following Garmin GPS devices:
 Etrex Summit, Etrex Vista Cx and GPSmap 60CSx. The function calls
 gpsbabel via the system. The gpsbabel program must be presen


Has anyone an explanation and a workaround ?


You will need to make the complete file available to us to diagnose 
this.  Is it in pgirmess 1.4.0?  Which topic?


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.


Re: [R] reference on contr.helmert and typo on its help page.

2009-11-08 Thread Peter Dalgaard

Gabor Grothendieck wrote:

On Sun, Nov 8, 2009 at 11:59 AM, Peng Yu  wrote:

On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch  wrote:

On 08/11/2009 11:03 AM, Peng Yu wrote:

I'm wondering which textbook discussed the various contrast matrices
mentioned in the help page of 'contr.helmert'. Could somebody let me
know?

Doesn't the reference on that page discuss them?

It does explain what the functions are. But I need a more basic and
complete reference. For example, I want to understand what 'Helmert
parametrization' (on page 33 of 'Statistical Models in S') is.



Just google for: Helmert contrasts


Or,

> contr.helmert(5)
  [,1] [,2] [,3] [,4]
1   -1   -1   -1   -1
21   -1   -1   -1
302   -1   -1
4003   -1
50004

> MASS::fractions(MASS::ginv(contr.helmert(5)))
 [,1]  [,2]  [,3]  [,4]  [,5]
[1,]  -1/2   1/2 0 0 0
[2,]  -1/6  -1/6   1/3 0 0
[3,] -1/12 -1/12 -1/12   1/4 0
[4,] -1/20 -1/20 -1/20 -1/20   1/5

and apply brains.

I.e., except for a slightly odd multiplier, the parameters represent the 
 difference between each level and the average of the preceding levels.


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

__
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] linear trend line and a quadratic trend line.

2009-11-08 Thread Eric Fail



Dear list users



How is it possible to visualise both a linear trend line and a quadratic trend 
line on a plot
of two variables?



Here my almost working exsample.



data(Duncan)

attach(Duncan)

plot(prestige ~ income)

abline(lm(prestige ~ income), col=2, lwd=2)



Now I would like to add yet another trend line, but this time a quadratic one. 
So I have two
trend lines. One linear trend line and a quadratic trend line. A trend line as 
if I had taken
I(income^2) into the equation.



I know I can make two models and compare them using anova, but for pedagogical 
resons I wold
like to visualise it. I know the trick from my past as an SPSS user, but I 
would like to do
this in R as well. Se it in SPSS
http://www.childrens-mercy.org/stats/weblog2006/QuadraticRegression.asp



Thanks in adcance



Eric

__
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] reference on contr.helmert and typo on its help page.

2009-11-08 Thread Gabor Grothendieck
On Sun, Nov 8, 2009 at 11:59 AM, Peng Yu  wrote:
> On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch  wrote:
>> On 08/11/2009 11:03 AM, Peng Yu wrote:
>>>
>>> I'm wondering which textbook discussed the various contrast matrices
>>> mentioned in the help page of 'contr.helmert'. Could somebody let me
>>> know?
>>
>> Doesn't the reference on that page discuss them?
>
> It does explain what the functions are. But I need a more basic and
> complete reference. For example, I want to understand what 'Helmert
> parametrization' (on page 33 of 'Statistical Models in S') is.
>

Just google for: Helmert contrasts

__
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] Counting non-empty levels of a factor

2009-11-08 Thread John Kane
With xx as your sample data will this work?  See ?addmargins

jj <- table(xx)

addmargins(jj, 2)

# or for both margins

addmargins(jj, c(1,2))

or 
apply(jj, 1, sum)


--- On Sun, 11/8/09, sylvain willart  wrote:

> From: sylvain willart 
> Subject: [R] Counting non-empty levels of a factor
> To: r-help@r-project.org
> Received: Sunday, November 8, 2009, 8:38 AM
> Hi everyone,
> 
> I'm struggling with a little problem for a while, and I'm
> wondering if
> anyone could help...
> 
> I have a dataset (from retailing industry) that indicates
> which brands
> are present in a panel of 500 stores,
> 
> store , brand
> 1 , B1
> 1 , B2
> 1 , B3
> 2 , B1
> 2 , B3
> 3 , B2
> 3 , B3
> 3 , B4
> 
> I would like to know how many brands are present in each
> store,
> 
> I tried:
> result <- aggregate(MyData$brand , by=list(MyData$store)
> , nlevels)
> 
> but I got:
> Group.1 x
> 1 , 4
> 2 , 4
> 3 , 4
> 
> which is not exactly the result I expected
> I would like to get sthg like:
> Group.1 x
> 1 , 3
> 2 , 2
> 3 , 3
> 
> Looking around, I found I can delete empty levels of factor
> using:
> problem.factor <- problem.factor[,drop=TRUE]
> But this solution isn't handy for me as I have many stores
> and should
> make a subset of my data for each store before dropping
> empty factor
> 
> I can't either counting the line for each store (N),
> because the same
> brand can appear several times in each store (several
> products for the
> same brand, and/or several weeks of observation)
> 
> I used to do this calculation using SAS with:
> proc freq data = MyData noprint ; by store ;
>  tables  brand / out = result ;
> run ;
> (the cool thing was I got a database I can merge with
> MyData)
> 
> any idea for doing that in R ?
> 
> Thanks in advance,
> 
> King Regards,
> 
> Sylvain Willart,
> PhD Marketing,
> IAE Lille, France
> 
> __
> 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.
> 


  __
Make your browsing faster, safer, and easier with the new Internet Explorer® 8. 
Opt
ernetexplorer/

__
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] rd doc truncated with R 2.10.0

2009-11-08 Thread Patrick Giraudoux

Hi,

I am routinely compiling a package and since I have moved to R 2.10.0, 
it troncates some section texts in the doc:


With the following section in the rd file:

\details{
The function calls gpsbabel via the system. The gpsbabel program must 
be present and on the user's PATH for the function to work see 
. The function has been tested on the 
following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and GPSmap 60CSx.

}

...compiling under R 2.9.2 (rcmd build --binary --auto-zip pgirmess) I 
get this


Details:

The function calls gpsbabel via the system. The gpsbabel program
must be present and on the user's PATH for the function to work,
see . The function has been tested on
the following Garmin GPS devices: Etrex Summit, Etrex Vista Cx and
GPSmap 60CSx.

and compiling now under R 2.10.0

Details:

The function has been tested on the following Garmin GPS devices:
Etrex Summit, Etrex Vista Cx and GPSmap 60CSx. The function calls
gpsbabel via the system. The gpsbabel program must be presen


Has anyone an explanation and a workaround ?

Best,

Patrick

__
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] reference on contr.helmert and typo on its help page.

2009-11-08 Thread Peng Yu
On Sun, Nov 8, 2009 at 10:11 AM, Duncan Murdoch  wrote:
> On 08/11/2009 11:03 AM, Peng Yu wrote:
>>
>> I'm wondering which textbook discussed the various contrast matrices
>> mentioned in the help page of 'contr.helmert'. Could somebody let me
>> know?
>
> Doesn't the reference on that page discuss them?

It does explain what the functions are. But I need a more basic and
complete reference. For example, I want to understand what 'Helmert
parametrization' (on page 33 of 'Statistical Models in S') is.

>> BTW, in R version 2.9.1, there is a typo on the help page of
>> 'contr.helmert' ('cont.helmert' should be 'contr.helmert').
>
> Thanks, I've fixed the typo.
>
> 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.


Re: [R] Strange Conflicts with Debian Repositories

2009-11-08 Thread Dirk Eddelbuettel
On Sun, Nov 08, 2009 at 10:31:12AM -0600, Dirk Eddelbuettel wrote:
> 
> On 8 November 2009 at 17:05, Lorenzo Isella wrote:
> | I am experiencing some really strange problems in updating my system
> | whenever I have both the R repository and the multimedia repository
> | available.
> [...]
> | deb http://cran.ch.r-project.org/bin/linux/debian lenny-cran/
> 
> This URLs doesn't exist / doesn't resolve. So don't use it. This one seems to
> work 
> 
>  http://stat.ethz.ch/CRAN/bin/linux/debian/lenny-cran/

Sorry, I meant

   http://stat.ethz.ch/CRAN/bin/linux/debian/  lenny-cran/

which works for me across the pond as well.

Dirk

-- 
Three out of two people have difficulties with fractions.

__
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] look up and Missing

2009-11-08 Thread John Kane
I'm not quite sure I understood the second queston but does this work?

subset(temp, xx$v2==-9)

subset(temp, xx$v2!= -9)

--- On Sun, 11/8/09, Ashta  wrote:

> From: Ashta 
> Subject: [R] look up and Missing
> To: r-h...@stat.math.ethz.ch
> Received: Sunday, November 8, 2009, 10:23 AM
> HI  R-Users
> 
> Assume that I have a data frame 'temp' with several
> variables (v1,v2,v3,v4,v5.).
> 
>   v1 v2 v3  v4 v5
>    1 
> 2   3   3    6
>    5  2  4    2 
>   0
>    2
> -9   5   4    3
>    6 
> 2   1   3    4
> 
> 1, I want to look at the entire row values of when v2 =-9
>    like
>          2
> -9   5   4    3
> 
> I wrote
> K<- list(if(temp$v2)==-9))
> 
> I wrote the like this but  it gave me  which is
> not correct.
>    False false false false false
> 
> 2. I want assign that values  as missing
> if   v2 = -9.  (ie., I want
> exclude from the analysis
> 
> How do I do it  in R?
> 
> Thanks in advance
> 
> __
> 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.
> 


  __
Yahoo! Canada Toolbar: Search from anywhere on the web, and bookmark your 
favourite sites. Download it now
http://ca.toolbar.yahoo.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] look up and Missing

2009-11-08 Thread David Winsemius


On Nov 8, 2009, at 11:08 AM, David Winsemius wrote:



On Nov 8, 2009, at 10:23 AM, Ashta wrote:


HI  R-Users

Assume that I have a data frame 'temp' with several variables  
(v1,v2,v3,v4,v5.).


v1 v2 v3  v4 v5
 1  2   3   36
 5  2  420
 2 -9   5   43
 6  2   1   34

1, I want to look at the entire row values of when v2 =-9
 like
   2 -9   5   43




I wrote
K<- list(if(temp$v2)==-9))


A further thought, that might be more useful if you were intending to  
supply a portion of a dataframe to an analytical function, would be  
the subset function:


t2 <- subset(temp, v2 != -9)

E. g.:

lm( v1 ~ v2 + v3, data= subset(temp, v2 != -9)




"if" would be the wrong R function to use. It's mostly for program  
control. And where did the "3" come from? You were working with the  
column temp$v2. Oh, you wanted a row rather than the column, "v2"?  
So how were you going to select that row? Perhaps:


K <-temp[ temp$v2 == -9, ]
K



I wrote the like this but  it gave me  which is not correct.
 False false false false false


I could not get your code to produce this. I got:
Error: unexpected '==' in "K<- list(if(temp$v2)=="



2. I want assign that values  as missing if   v2 = -9.  (ie., I want
exclude from the analysis

How do I do it  in R?


Your request is not well specified at least to my reading, because I  
could not tell if you wanted the re-assignment to occur in temp (and  
that was after I came down on the row side of the whether you wanted  
a row or column.) . The following assumes you wanted the row in  
question (created above)  modified outside of "temp".


> is.na(K) <- K == -9
> K
 v1 v2 v3 v4 v5
3  2 NA  5  4  3

If you had used ifelse you would have gotten close, but the data  
type would have been a list, which may not have been what you  
expected:


> K <- ifelse(K==-9, NA, K)
> K
[[1]]
[1] 2

[[2]]
[1] NA

[[3]]
[1] 5

[[4]]
[1] 4

[[5]]
[1] 3




--

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] look up and Missing

2009-11-08 Thread Jorge Ivan Velez
Dear Ashta,

Is this what you want?

x <- read.table(textConnection("v1 v2 v3  v4 v5
  1  2   3   36
  5  2  420
  2 -9   5   43
  6  2   1   34"), header = TRUE)
closeAllConnections()
x

# Option 1
x1 <- x  # copy of x just for this example
x1$v2[which(x1$v2 == -9)] <- NA
x1

# Option 2
x2 <- x  # copy of x just for this example
x2$v2 <- with(x2, ifelse(v2 == -9, NA, v2))
x2

HTH,
Jorge



On Sun, Nov 8, 2009 at 10:23 AM, Ashta <> wrote:

> HI  R-Users
>
> Assume that I have a data frame 'temp' with several variables
> (v1,v2,v3,v4,v5.).
>
>  v1 v2 v3  v4 v5
>   1  2   3   36
>   5  2  420
>   2 -9   5   43
>   6  2   1   34
>
> 1, I want to look at the entire row values of when v2 =-9
>   like
> 2 -9   5   43
>
> I wrote
> K<- list(if(temp$v2)==-9))
>
> I wrote the like this but  it gave me  which is not correct.
>   False false false false false
>
> 2. I want assign that values  as missing if   v2 = -9.  (ie., I want
> exclude from the analysis
>
> How do I do it  in R?
>
> Thanks in advance
>
> __
> 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] look up and Missing

2009-11-08 Thread Dimitris Rizopoulos

try this:

temp.new <- temp[temp$v2 != -9, ]
temp.new


I hope it helps.

Best,
Dimitris


Ashta wrote:

HI  R-Users

Assume that I have a data frame 'temp' with several variables (v1,v2,v3,v4,v5.).

  v1 v2 v3  v4 v5
   1  2   3   36
   5  2  420
   2 -9   5   43
   6  2   1   34

1, I want to look at the entire row values of when v2 =-9
   like
 2 -9   5   43

I wrote
K<- list(if(temp$v2)==-9))

I wrote the like this but  it gave me  which is not correct.
   False false false false false

2. I want assign that values  as missing if   v2 = -9.  (ie., I want
exclude from the analysis

How do I do it  in R?

Thanks in advance

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



--
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014

__
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] look up and Missing

2009-11-08 Thread David Winsemius


On Nov 8, 2009, at 10:23 AM, Ashta wrote:


HI  R-Users

Assume that I have a data frame 'temp' with several variables  
(v1,v2,v3,v4,v5.).


 v1 v2 v3  v4 v5
  1  2   3   36
  5  2  420
  2 -9   5   43
  6  2   1   34

1, I want to look at the entire row values of when v2 =-9
  like
2 -9   5   43




I wrote
K<- list(if(temp$v2)==-9))


"if" would be the wrong R function to use. It's mostly for program  
control. And where did the "3" come from? You were working with the  
column temp$v2. Oh, you wanted a row rather than the column, "v2"? So  
how were you going to select that row? Perhaps:


K <-temp[ temp$v2 == -9, ]
K



I wrote the like this but  it gave me  which is not correct.
  False false false false false


I could not get your code to produce this. I got:
Error: unexpected '==' in "K<- list(if(temp$v2)=="



2. I want assign that values  as missing if   v2 = -9.  (ie., I want
exclude from the analysis

How do I do it  in R?


Your request is not well specified at least to my reading, because I  
could not tell if you wanted the re-assignment to occur in temp (and  
that was after I came down on the row side of the whether you wanted a  
row or column.) . The following assumes you wanted the row in question  
(created above)  modified outside of "temp".


> is.na(K) <- K == -9
> K
  v1 v2 v3 v4 v5
3  2 NA  5  4  3

If you had used ifelse you would have gotten close, but the data type  
would have been a list, which may not have been what you expected:


> K <- ifelse(K==-9, NA, K)
> K
[[1]]
[1] 2

[[2]]
[1] NA

[[3]]
[1] 5

[[4]]
[1] 4

[[5]]
[1] 3




--

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] Strange Conflicts with Debian Repositories

2009-11-08 Thread Dirk Eddelbuettel

On 8 November 2009 at 17:05, Lorenzo Isella wrote:
| I am experiencing some really strange problems in updating my system
| whenever I have both the R repository and the multimedia repository
| available.
[...]
| deb http://cran.ch.r-project.org/bin/linux/debian lenny-cran/

This URLs doesn't exist / doesn't resolve. So don't use it. This one seems to
work 

 http://stat.ethz.ch/CRAN/bin/linux/debian/lenny-cran/

so try that instead. Likewise for the sources. Also note the trailing slash.

This question, by the way, would be much more appropriate for the
r-sig-debian list which is dedicated to Debian / Ubuntu questions.

Dirk

-- 
Three out of two people have difficulties with fractions.

__
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] reference on contr.helmert and typo on its help page.

2009-11-08 Thread David Winsemius


On Nov 8, 2009, at 11:03 AM, Peng Yu wrote:


I'm wondering which textbook discussed the various contrast matrices
mentioned in the help page of 'contr.helmert'. Could somebody let me
know?



My version of "Modern Applied Statistics in S" (aka MASS) deals with  
it in enough detail for a person with background to understand. At one  
point I asked Ripley whether later editions of MASS expanded the  
coverage of that topic but my memory is that he did not reply. Coming  
from a perspective that emphasized regression, I found it rather terse  
(2 pages), so there might be one of the more recently published texts  
that spends more time developing the concepts from basics.


"Statistical Models in S" also covers the topic (5 pages)  in an early  
chapter, again at a level that assumes prior training in ANOVA models.



BTW, in R version 2.9.1, there is a typo on the help page of
'contr.helmert' ('cont.helmert' should be 'contr.helmert').

__

--

David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] look up and Missing

2009-11-08 Thread jim holtman
Here is how to find out which rows contain -9 and then you can do with
it as you please:

> x
  v1 v2 v3 v4 v5
1  1  2  3  3  6
2  5  2  4  2  0
3  2 -9  5  4  3
4  6  2  1  3  4
> which(apply(x, 1, function(.row) any(.row == -9)))
[1] 3
>
>


On Sun, Nov 8, 2009 at 10:23 AM, Ashta  wrote:
> HI  R-Users
>
> Assume that I have a data frame 'temp' with several variables 
> (v1,v2,v3,v4,v5.).
>
>  v1 v2 v3  v4 v5
>   1  2   3   3    6
>   5  2  4    2    0
>   2 -9   5   4    3
>   6  2   1   3    4
>
> 1, I want to look at the entire row values of when v2 =-9
>   like
>         2 -9   5   4    3
>
> I wrote
> K<- list(if(temp$v2)==-9))
>
> I wrote the like this but  it gave me  which is not correct.
>   False false false false false
>
> 2. I want assign that values  as missing if   v2 = -9.  (ie., I want
> exclude from the analysis
>
> How do I do it  in R?
>
> Thanks in advance
>
> __
> 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 that 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] Recall: Scheffe test

2009-11-08 Thread Mangalani Peter Makananisa
Mangalani Peter Makananisa would like to recall the message, "Scheffe test".

Please Note: This email and its contents are subject to our email legal notice 
which can be viewed at http://www.sars.gov.za/Email_Disclaimer.pdf 

[[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] Scheffe test

2009-11-08 Thread Mangalani Peter Makananisa
 

Dear Prof,

 

Please help me with the R code which compute SCHEFFE TEST

 

Thanking you in advance

 

Kind regards

 

Mangalani Peter Makananisa

Statistical Analyst

South African Revenue Service (SARS)

Segmentation and Research : Data Modelling

Tel: +2712 422 7357

Cell: +2782 456 4669

Fax:+2712 422 6579

 

 


Please Note: This email and its contents are subject to our email legal notice 
which can be viewed at http://www.sars.gov.za/Email_Disclaimer.pdf 

[[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] Scheffe test

2009-11-08 Thread Mangalani Peter Makananisa
 

Dear all,

 

Please help me with the R code which compute SCHEFFE TEST

 

Thanking you in advance

 

Kind regards

 

Mangalani Peter Makananisa

Statistical Analyst

South African Revenue Service (SARS)

Segmentation and Research : Data Modelling

Tel: +2712 422 7357

Cell: +2782 456 4669

Fax:+2712 422 6579

 

 


Please Note: This email and its contents are subject to our email legal notice 
which can be viewed at http://www.sars.gov.za/Email_Disclaimer.pdf 

[[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] Scheffe test

2009-11-08 Thread Mangalani Peter Makananisa
 

Dear all,

 

Please help me with the R code which compute SCHEFFE TEST

 

Thanking you in advance

 

Kind regards

 

Mangalani Peter Makananisa

Statistical Analyst

South African Revenue Service (SARS)

Segmentation and Research : Data Modelling

Tel: +2712 422 7357

Cell: +2782 456 4669

Fax:+2712 422 6579

 

 


Please Note: This email and its contents are subject to our email legal notice 
which can be viewed at http://www.sars.gov.za/Email_Disclaimer.pdf 
__
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] Strange Conflicts with Debian Repositories

2009-11-08 Thread Lorenzo Isella
Dear All,
I have just installed a fresh Debian testing (squeeze) on my system
(amd64 architecture).
I am experiencing some really strange problems in updating my system
whenever I have both the R repository and the multimedia repository
available.
This is my source.list (when I disable the multimedia repository)

~$ cat /etc/apt/sources.list
deb http://ftp.ch.debian.org/debian/ testing main contrib non-free
deb-src http://ftp.ch.debian.org/debian/ testing main non-free contrib

deb http://security.debian.org/ testing/updates main contrib non-free
deb-src http://security.debian.org/ testing/updates main contrib non-free

#deb http://www.debian-multimedia.org testing main

#deb http:///bin/linux/debian lenny-cran/
#deb-src http:///bin/linux/debian lenny-cran/

# Add R-repository
deb http://cran.ch.r-project.org/bin/linux/debian lenny-cran/
deb-src http://cran.ch.r-project.org/debian/bin/linux/debian lenny-cran/

This way I have no trouble to update/upgrade my system (only, I do not
see the sources of the R packages, but this is really a minor issue
for now)

~$ sudo apt-get update
Get:1 http://security.debian.org testing/updates Release.gpg [835B]
Ign http://security.debian.org testing/updates/main Translation-en_US
Ign http://security.debian.org testing/updates/contrib Translation-en_US
Get:2 http://ftp.ch.debian.org testing Release.gpg [835B]
Ign http://ftp.ch.debian.org testing/main Translation-en_US
Ign http://ftp.ch.debian.org testing/contrib Translation-en_US
Ign http://security.debian.org testing/updates/non-free Translation-en_US
Hit http://security.debian.org testing/updates Release
Ign http://ftp.ch.debian.org testing/non-free Translation-en_US
Hit http://ftp.ch.debian.org testing Release
Ign http://security.debian.org testing/updates/main Packages/DiffIndex
Hit http://ftp.ch.debian.org testing/main Packages/DiffIndex
Ign http://security.debian.org testing/updates/contrib Packages/DiffIndex
Ign http://security.debian.org testing/updates/non-free Packages/DiffIndex
Ign http://security.debian.org testing/updates/main Sources/DiffIndex
Ign http://security.debian.org testing/updates/contrib Sources/DiffIndex
Ign http://security.debian.org testing/updates/non-free Sources/DiffIndex
Get:3 http://cran.ch.r-project.org lenny-cran/ Release.gpg [197B]
Ign http://cran.ch.r-project.org lenny-cran/ Translation-en_US
Ign http://cran.ch.r-project.org lenny-cran/ Release.gpg
Hit http://ftp.ch.debian.org testing/contrib Packages/DiffIndex
Ign http://ftp.ch.debian.org testing/non-free Packages/DiffIndex
Hit http://ftp.ch.debian.org testing/main Sources/DiffIndex
Ign http://ftp.ch.debian.org testing/non-free Sources/DiffIndex
Hit http://ftp.ch.debian.org testing/contrib Sources/DiffIndex
Hit http://security.debian.org testing/updates/main Packages
Hit http://security.debian.org testing/updates/contrib Packages
Hit http://security.debian.org testing/updates/non-free Packages
Hit http://security.debian.org testing/updates/main Sources
Hit http://security.debian.org testing/updates/contrib Sources
Get:4 http://cran.ch.r-project.org lenny-cran/ Release [1,288B]
Hit http://ftp.ch.debian.org testing/non-free Packages
Hit http://ftp.ch.debian.org testing/non-free Sources
Hit http://security.debian.org testing/updates/non-free Sources
Ign http://cran.ch.r-project.org lenny-cran/ Release
Ign http://cran.ch.r-project.org lenny-cran/ Packages
Ign http://cran.ch.r-project.org lenny-cran/ Sources
Ign http://cran.ch.r-project.org lenny-cran/ Packages
Ign http://cran.ch.r-project.org lenny-cran/ Sources
Get:5 http://cran.ch.r-project.org lenny-cran/ Packages [11.4kB]
Err http://cran.ch.r-project.org lenny-cran/ Sources
  404  Not Found
Fetched 14.5kB in 3s (3,851B/s)
W: Failed to fetch
http://cran.ch.r-project.org/debian/bin/linux/debian/lenny-cran/Sources.gz
 404  Not Found

E: Some index files failed to download, they have been ignored, or old
ones used instead.

However (and this is really odd...) if I uncomment the multimedia
repository this is what happens

$ sudo apt-get update
Err http://www.debian-multimedia.org testing Release.gpg
  Could not resolve 'www.debian-multimedia.org'
Err http://security.debian.org testing/updates Release.gpg
  Could not resolve 'security.debian.org'
Err http://security.debian.org testing/updates/main Translation-en_US
  Could not resolve 'security.debian.org'
Err http://www.debian-multimedia.org testing/main Translation-en_US
  Could not resolve 'www.debian-multimedia.org'
Err http://security.debian.org testing/updates/contrib
Translation-en_US
  Could not resolve 'security.debian.org'
Err http://security.debian.org testing/updates/non-free Translation-en_US
  Could not resolve 'security.debian.org'
Err http://cran.ch.r-project.org lenny-cran/ Release.gpg
  Could not resolve 'cran.ch.r-project.org'
Err http://cran.ch.r-project.org lenny-cran/ Translation-en_US
  Could not resolve 'cran.ch.r-project.org'
Err http://cran.ch.r-project.org lenny-cran/ Release.gpg
  Could not resolve 'cran.ch.r-proje

[R] reference on contr.helmert and typo on its help page.

2009-11-08 Thread Peng Yu
I'm wondering which textbook discussed the various contrast matrices
mentioned in the help page of 'contr.helmert'. Could somebody let me
know?

BTW, in R version 2.9.1, there is a typo on the help page of
'contr.helmert' ('cont.helmert' should be 'contr.helmert').

__
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] Counting non-empty levels of a factor

2009-11-08 Thread sylvain willart
Thanks a lot for those solutions,
Both are working great, and they do slightly different (but both very
interesting) things,
Moreover, I learned about the length() function ... one more to add to
my personal cheat sheet
King Regards

2009/11/8 David Winsemius :
>
> On Nov 8, 2009, at 9:11 AM, David Winsemius wrote:
>
>>
>> On Nov 8, 2009, at 8:38 AM, sylvain willart wrote:
>>
>>> Hi everyone,h
>>>
>>> I'm struggling with a little problem for a while, and I'm wondering if
>>> anyone could help...
>>>
>>> I have a dataset (from retailing industry) that indicates which brands
>>> are present in a panel of 500 stores,
>>>
>>> store , brand
>>> 1 , B1
>>> 1 , B2
>>> 1 , B3
>>> 2 , B1
>>> 2 , B3
>>> 3 , B2
>>> 3 , B3
>>> 3 , B4
>>>
>>> I would like to know how many brands are present in each store,
>>>
>>> I tried:
>>> result <- aggregate(MyData$brand , by=list(MyData$store) , nlevels)
>>>
>>> but I got:
>>> Group.1 x
>>> 1 , 4
>>> 2 , 4
>>> 3 , 4
>>>
>>> which is not exactly the result I expected
>>> I would like to get sthg like:
>>> Group.1 x
>>> 1 , 3
>>> 2 , 2
>>> 3 , 3
>>
>> Try:
>>
>> result <- aggregate(MyData$brand , by=list(MyData$store) , length)
>>
>> Quick, easy and generalizes to other situations. The factor levels got
>> carried along identically, but length counts the number of elements in the
>> list returned by tapply.
>
> Which may not have been what you asked for as this would demonstrate. You
> probably wnat the second solution:
> mydata2 <- rbind(MyData, MyData)
>> result <- aggregate(mydata2$brand , by=list(mydata2$store) , length)
>> result
>  Group.1 x
> 1       1 6
> 2       2 4
> 3       3 6
>
>> result <- aggregate(mydata2$brand , by=list(mydata2$store) , function(x)
>> nlevels(factor(x)))
>> result
>  Group.1 x
> 1       1 3
> 2       2 2
> 3       3 3
>
>>>
>>> Looking around, I found I can delete empty levels of factor using:
>>> problem.factor <- problem.factor[,drop=TRUE]
>>
>> If you reapply the function, factor, you get the same result. So you could
>> have done this:
>>
>> > result <- aggregate(MyData$brand , by=list(MyData$store) , function(x)
>> > nlevels(factor(x)))
>> > result
>>  Group.1 x
>> 1       1 3
>> 2       2 2
>> 3       3 3
>>
>>
>>
>>> But this solution isn't handy for me as I have many stores and should
>>> make a subset of my data for each store before dropping empty factor
>>>
>>> I can't either counting the line for each store (N), because the same
>>> brand can appear several times in each store (several products for the
>>> same brand, and/or several weeks of observation)
>>>
>>> I used to do this calculation using SAS with:
>>> proc freq data = MyData noprint ; by store ;
>>> tables  brand / out = result ;
>>> run ;
>>> (the cool thing was I got a database I can merge with MyData)
>>>
>>> any idea for doing that in R ?
>>>
>>> Thanks in advance,
>>>
>>> King Regards,
>>>
>>> Sylvain Willart,
>>> PhD Marketing,
>>> IAE Lille, France
>>>
>>> __
>>> 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.
>>
>> David Winsemius, MD
>> Heritage Laboratories
>> West Hartford, CT
>>
>> __
>> 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.
>
> David Winsemius, MD
> Heritage Laboratories
> West Hartford, CT
>
>

__
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] save an object by dynamicly created name

2009-11-08 Thread Hao Cen
Hi Henrik,

I am using your saveObject/loadObject to handle over 1000 matrices. It
worked beautifully. Because I need to load those matrices often for
evaluating a few functions on them and those matrices do not fit all in
memory at once, is there a way to speed up the loading part? I tried save
all the binary files to /dev/shm  (shared memory section in linux) but the
speed of loadObject on /dev/shm remains the same as on the disk.

Thanks

Hao



-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On
Behalf Of Henrik Bengtsson
Sent: Monday, November 02, 2009 12:34 AM
To: David Winsemius
Cc: r-help@r-project.org; jeffc
Subject: Re: [R] save an object by dynamicly created name

On Sun, Nov 1, 2009 at 9:18 PM, David Winsemius 
wrote:
>
> On Nov 1, 2009, at 11:28 PM, Henrik Bengtsson wrote:
>
>> On Sun, Nov 1, 2009 at 7:48 PM, David Winsemius 
>> wrote:
>>>
>>> On Nov 1, 2009, at 10:16 PM, Henrik Bengtsson wrote:
>>>
 path <- "data";
 dir.create(path);

 for (i in 1:10) {
  m <- i:5;
  filename <- sprintf("m%02d.Rbin", i);
  pathname <- file.path(path, filename);
  save(m, file=pathname);
 }

>>>
>>> That would result in each of the ten files containing an object with the
>>> same  name == "m". (Also on my system R data files have type Rdta.) So I
>>> thought what was requested might have been a slight mod:
>>>
>>> path <- "~/";
>>> dir.create(path);
>>>
>>> for (i in 1:10) {
>>>  assign( paste("m", i, sep=""),  i:5)
>>>  filename <- sprintf("m%02d.Rdta", i)
>>>  pathname <- file.path(path, filename)
>>>  obj =get(paste("m", i, sep=""))
>>>  save(obj, file=pathname)
>>> }
>>
>> Then a more convenient solution is to use saveObject() and
>> loadObject() of R.utils.  saveObject() does not save the name of the
>> object save.
>
> The OP asked for this outcome :
>
> " I would like to save m as m1, m2, m3 ...,
> to file /home/data/m1, /home/data/m2, home/data/m3, ..."
>
>
>>  If you want to save multiple objects, the wrap them up
>> in a list.
>
> I agree that a list would makes sense if it were to be stored in one file
,
> although it was not what requested.

That comment was not for the OP, but for saveObject()/loadObject() in
general.

> But wouldn't that require assign()-ing a name before list()-wrapping?

Nope, the whole point of using saveObject()/loadObject() is to save
the objects/values without their names that you happens to choose in
the current session, and to avoid overwriting existing ones in your
next session. My example could also have been:

library("R.utils");
saveObject(list(a=1,b=LETTERS,c=Sys.time()), file="foo.Rbin");
y <- loadObject("foo.Rbin");
z <- loadObject("foo.Rbin");
stopifnot(identical(y,z));

If you really want to attach the elements of the saved list, do:

attachLocally(loadObject("foo.Rbin"));
> str(a)
 num 1
> str(b)
 chr [1:26] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" ...
> str(c)
 POSIXct[1:1], format: "2009-11-01 21:30:41"

>
> I suppose we ought to mention that the use of assign to create a variable
is
> a FAQ ... 7.21? Yep, I have now referred to it a sufficient number of
times
> to refer to it by number.
>
>
http://cran.r-project.org/doc/FAQ/R-FAQ.html#How-can-I-turn-a-string-into-a-
variable_003f

My personal take on assign() and get() is that if you find yourself
using them (at this level), there is a good chance there exists a
better solution that you should use instead.

My $.02

/H

>
> --
> David
>
>>  loadObject() does not assign variable, but instead return
>> them. Example:
>>
>> library("R.utils");
>> x <- list(a=1,b=LETTERS,c=Sys.time());
>> saveObject(x, file="foo.Rbin");
>> y <- loadObject("foo.Rbin");
>> stopifnot(identical(x,y));
>
>>
>> So, for the original example, I'd recommend:
>>
>> library("R.utils");
>> path <- "data";
>> mkdirs(path);
>>
>> for (i in 1:10) {
>>  m <- i:5;
>>  filename <- sprintf("m%02d.Rbin", i);
>>  saveObject(m, file=filename, path=path);
>> }
>>
>> and loading the objects back as:
>>
>> for (i in 1:10) {
>>  filename <- sprintf("m%02d.Rbin", i);
>>  m <- loadObject(filename, path=path);
>>  print(m);
>> }
>> /Henrik
>>
>>>
>>> --
>>> David.
>>>
 /H

 On Sun, Nov 1, 2009 at 6:53 PM, jeffc  wrote:
>
> Hi,
>
> I would like to save a few dynamically created objects to disk. The
> following is the basic flow of the code segment
>
> for(i = 1:10) {
>  m = i:5
>  save(m, file = ...) ## ???
> }
> To distinguish different objects to be saved, I would like to save m
as
> m1,
> m2, m3 ..., to file /home/data/m1, /home/data/m2, home/data/m3, ...
>
> I tried a couple of methods on translating between object names and
> strings
> (below) but couldn't get it to work.
> https://stat.ethz.ch/pipermail/r-help/2008-November/178965.html
> http://tolstoy.newcastle.edu.au/R/help/04/08/2673.html
>
> Any suggestions would be appreciated.
>
> thanks
>
> Hao
>
>

[R] how to remove one of any two indices with a correlation greater than 0.90 in a matrix (or data.frame)

2009-11-08 Thread bbslover

my code is not right below:
rm(list=ls())
#define data.frame
a=c(1,2,3,5,6); b=c(1,2,3,4,7); c=c(1,2,3,4,8); d=c(1,2,3,5,1);
e=c(1,2,3,5,7) 
data.f=data.frame(a,b,c,d,e) 
#backup data.f
origin.data<-data.f 
#get correlation matrix
cor.matrix<-cor(origin.data) 
#backup correlation matrix
origin.cor<-cor.matrix
#get dim
dim.cor<-dim(origin.cor) 
#perform Loop
n<-0
for(i in 1:(dim.cor[1]-1)) 
{ 
  for(j in (i+1):(dim.cor[2]))
   { 
  if (cor.matrix[i,j]>=0.95) 
  { 
  data.f<-data.f[,-(i-n)] 
  n<-1
  break
  } 
} 
} 
origin.cor 
origin.data
data.f 
cor(data.f)

how write the code to do with my questions? and have a simple way?  
-- 
View this message in context: 
http://old.nabble.com/how-to-remove-one-of-any-two-indices-with-a-correlation-greater-than-0.90-in-a-matrix-%28or-data.frame%29-tp26254082p26254082.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] Windows 7 editor - I can't make RWinEdt work

2009-11-08 Thread Uwe Ligges

Aha, what is that blog post and what does not work  for you?
I haven't got any report so far and do not have Windows 7 easily 
available yet.


Best,
Uwe Ligges



Peter Flom wrote:

Good morning

I just got a new computer with Windows 7.  R works fine, but the editor I am used to 
using "RWinEdt" does not.  I did find one blog post on how to get RWinEdt to 
work in Windows 7, but I could not get those instructions to work either.

Is there a patch for RWinEdt?

If not, is there another good R editor that works under Windows 7?

I tried RSiteSearch with various combinations of Windows 7 and Editor and so 
on, but found nothing.  I also tried googling on these terms.

Thanks

Peter

Peter L. Flom, PhD
Statistical Consultant
Website: www DOT peterflomconsulting DOT com
Writing; http://www.associatedcontent.com/user/582880/peter_flom.html
Twitter:   @peterflom

__
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] look up and Missing

2009-11-08 Thread Ashta
HI  R-Users

Assume that I have a data frame 'temp' with several variables (v1,v2,v3,v4,v5.).

  v1 v2 v3  v4 v5
   1  2   3   36
   5  2  420
   2 -9   5   43
   6  2   1   34

1, I want to look at the entire row values of when v2 =-9
   like
 2 -9   5   43

I wrote
K<- list(if(temp$v2)==-9))

I wrote the like this but  it gave me  which is not correct.
   False false false false false

2. I want assign that values  as missing if   v2 = -9.  (ie., I want
exclude from the analysis

How do I do it  in R?

Thanks in advance

__
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] Windows 7 editor - I can't make RWinEdt work

2009-11-08 Thread Peter Flom
Good morning

I just got a new computer with Windows 7.  R works fine, but the editor I am 
used to using "RWinEdt" does not.  I did find one blog post on how to get 
RWinEdt to work in Windows 7, but I could not get those instructions to work 
either.

Is there a patch for RWinEdt?

If not, is there another good R editor that works under Windows 7?

I tried RSiteSearch with various combinations of Windows 7 and Editor and so 
on, but found nothing.  I also tried googling on these terms.

Thanks

Peter

Peter L. Flom, PhD
Statistical Consultant
Website: www DOT peterflomconsulting DOT com
Writing; http://www.associatedcontent.com/user/582880/peter_flom.html
Twitter:   @peterflom

__
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] Counting non-empty levels of a factor

2009-11-08 Thread David Winsemius


On Nov 8, 2009, at 9:11 AM, David Winsemius wrote:



On Nov 8, 2009, at 8:38 AM, sylvain willart wrote:


Hi everyone,

I'm struggling with a little problem for a while, and I'm wondering  
if

anyone could help...

I have a dataset (from retailing industry) that indicates which  
brands

are present in a panel of 500 stores,

store , brand
1 , B1
1 , B2
1 , B3
2 , B1
2 , B3
3 , B2
3 , B3
3 , B4

I would like to know how many brands are present in each store,

I tried:
result <- aggregate(MyData$brand , by=list(MyData$store) , nlevels)

but I got:
Group.1 x
1 , 4
2 , 4
3 , 4

which is not exactly the result I expected
I would like to get sthg like:
Group.1 x
1 , 3
2 , 2
3 , 3


Try:

result <- aggregate(MyData$brand , by=list(MyData$store) , length)

Quick, easy and generalizes to other situations. The factor levels  
got carried along identically, but length counts the number of  
elements in the list returned by tapply.


Which may not have been what you asked for as this would demonstrate.  
You probably wnat the second solution:

mydata2 <- rbind(MyData, MyData)
> result <- aggregate(mydata2$brand , by=list(mydata2$store) , length)
> result
  Group.1 x
1   1 6
2   2 4
3   3 6

> result <- aggregate(mydata2$brand , by=list(mydata2$store) ,  
function(x) nlevels(factor(x)))

> result
  Group.1 x
1   1 3
2   2 2
3   3 3



Looking around, I found I can delete empty levels of factor using:
problem.factor <- problem.factor[,drop=TRUE]


If you reapply the function, factor, you get the same result. So you  
could have done this:


> result <- aggregate(MyData$brand , by=list(MyData$store) ,  
function(x) nlevels(factor(x)))

> result
 Group.1 x
1   1 3
2   2 2
3   3 3




But this solution isn't handy for me as I have many stores and should
make a subset of my data for each store before dropping empty factor

I can't either counting the line for each store (N), because the same
brand can appear several times in each store (several products for  
the

same brand, and/or several weeks of observation)

I used to do this calculation using SAS with:
proc freq data = MyData noprint ; by store ;
tables  brand / out = result ;
run ;
(the cool thing was I got a database I can merge with MyData)

any idea for doing that in R ?

Thanks in advance,

King Regards,

Sylvain Willart,
PhD Marketing,
IAE Lille, France

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] plot.window arguments being ignored?

2009-11-08 Thread Bryan Hanson
Dan, with base graphics, the> plot command determines the ylim, >lines can
only add lines/points to an exisiting graph, it cannot modify the existing
ylim.

You have two choices.  1) Before proceeding to the graph, determine which of
the two data sets has the greater range and call >plot with that range, then
the subsequent >lines command will "fit".  2) use a dlfferent graphics
system such as lattice or ggplot2.  In these packages, the entire graph is
computed in the background, then plotted explicitly.  So as you build up the
graph by adding things to it, the routines automatically resize as
necessary.  Both lattice and ggplot2 have books and very nice web sites
where you can find an example like yours.

Bryan
*
Bryan Hanson
Acting Chair
Professor of Chemistry & Biochemistry
DePauw University, Greencastle IN USA



On 11/7/09 6:28 PM, "AR_user"  wrote:

> 
> Hi,
> Why, when I run the script below, is my y-axis range command being ignored?
> I.e., the y axis graphed consistently fits the line specified in the plot
> command, but never fits the line I'm trying to add in the subsequent line
> command.
> This is my first R script, so if I'm doing anything else wacky that jumps
> out, please let me know.  I'm an advanced Stata user, but R is new to me.
> The script is designed to simulate a binary-choice decision field theory
> example with internally-controlled stop time, as per Neural Networks 19
> (2006) 1047-1058.
> Thanks in advance!!
> -Dan
> 
> 
> #-
> # Starting preference state; two choices
> preference<-c(0,0)
> # Preference threshold; when one of the preference states exceeds this
> threshhold, deliberation stops
> theta<- 2
> # Contrast matrix
> contrast<-array(c(1, -1, -1, 1), dim=c(2,2))
> # Value Matrix; three possible states of the world
> value<-array(c(1,0, 0, 1, .5,.5), dim=c(2, 3))
> # Feedback Matrix
> feedback<-array(c(1, .5, .5, 1), dim=c(2,2))
> 
> # Time
> t<- array(0, dim=c(1,1))
> 
> # Arrays to keep track of history for graphing
> time_for_graph<-array(t, dim=c(1,1))
> preference_for_graph<-array(preference, dim=c(1,2))
> 
> # Moving through time until preference crosses the threshold theta
> while (preference[1] 
> #stochastic weights for this time period
> weight<-rnorm(3,0,1)
> 
> #updating preference
> preference<-feedback%*%preference + contrast%*%value%*%weight
> 
> #updating history
> t<-t+1
> time_for_graph<-rbind(time_for_graph, t)
> preference_for_graph<-rbind(preference_for_graph, array(preference,
> dim=c(1,2)))
> 
> #updating graph ranges
> ry<-range(preference_for_graph)
> rx<-range(time_for_graph)
> 
> #updating graph 
> plot(time_for_graph[,1], preference_for_graph[,1], type="b",
> plot.window(rx, ry), xlab="Time", ylab="Preference")
> lines(time_for_graph[,1], preference_for_graph[,2], type="b", col="red",
> ylim=ry)
> }
> #-

__
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] Counting non-empty levels of a factor

2009-11-08 Thread David Winsemius


On Nov 8, 2009, at 8:38 AM, sylvain willart wrote:


Hi everyone,

I'm struggling with a little problem for a while, and I'm wondering if
anyone could help...

I have a dataset (from retailing industry) that indicates which brands
are present in a panel of 500 stores,

store , brand
1 , B1
1 , B2
1 , B3
2 , B1
2 , B3
3 , B2
3 , B3
3 , B4

I would like to know how many brands are present in each store,

I tried:
result <- aggregate(MyData$brand , by=list(MyData$store) , nlevels)

but I got:
Group.1 x
1 , 4
2 , 4
3 , 4

which is not exactly the result I expected
I would like to get sthg like:
Group.1 x
1 , 3
2 , 2
3 , 3


Try:

result <- aggregate(MyData$brand , by=list(MyData$store) , length)

Quick, easy and generalizes to other situations. The factor levels got  
carried along identically, but length counts the number of elements in  
the list returned by tapply.


Looking around, I found I can delete empty levels of factor using:
problem.factor <- problem.factor[,drop=TRUE]


If you reapply the function, factor, you get the same result. So you  
could have done this:


> result <- aggregate(MyData$brand , by=list(MyData$store) ,  
function(x) nlevels(factor(x)))

> result
  Group.1 x
1   1 3
2   2 2
3   3 3




But this solution isn't handy for me as I have many stores and should
make a subset of my data for each store before dropping empty factor

I can't either counting the line for each store (N), because the same
brand can appear several times in each store (several products for the
same brand, and/or several weeks of observation)

I used to do this calculation using SAS with:
proc freq data = MyData noprint ; by store ;
tables  brand / out = result ;
run ;
(the cool thing was I got a database I can merge with MyData)

any idea for doing that in R ?

Thanks in advance,

King Regards,

Sylvain Willart,
PhD Marketing,
IAE Lille, France

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


David Winsemius, MD
Heritage Laboratories
West Hartford, CT

__
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] influence.measures(stats): hatvalues(model, ...)

2009-11-08 Thread Viechtbauer Wolfgang (STAT)
Not sure what you mean.

yi <- c(2,3,2,4,3,6)
xi <- c(1,4,3,2,4,5)

res <- lm(yi ~ xi)
hatvalues(res)

X <- cbind(1, xi)
diag( X%*%solve(t(X)%*%X)%*%t(X) )

Same result.

Best,

--
Wolfgang Viechtbauerhttp://www.wvbauer.com/
Department of Methodology and StatisticsTel: +31 (0)43 388-2277
School for Public Health and Primary Care   Office Location:
Maastricht University, P.O. Box 616 Room B2.01 (second floor)
6200 MD Maastricht, The Netherlands Debyeplein 1 (Randwyck)

From: r-help-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf Of 
Sigmund Freud [ss_freud...@yahoo.com]
Sent: Sunday, November 08, 2009 8:14 AM
To: r-help@r-project.org
Subject: [R]  influence.measures(stats): hatvalues(model, ...)

Hello:

I am trying to understand the method 'hatvalues(...)', which returns something 
similar to the diagonals of the plain vanilla hat matrix [X(X'X)^(-1)X'], but 
not quite.

A Fortran programmer I am not, but tracing through the code it looks like 
perhaps some sort of correction based on the notion of 'leave-one-out' variance 
is being applied.

Whatever the difference, in simulations 'hatvalues' appears to perform much 
better in the context of identifying outiers using Cook's Distance than the 
diagonals of the plain vanilla hat matrix. (As in 
http://en.wikipedia.org/wiki/Cook's_distance).

Would prefer to understand a little more when using this method. I have 
downloaded the freely available references cited in the help and am in the 
process of digesting them. If someone with knowledge could offer a pointer on 
the most efficient way to get at why 'hatvalues' does what it does, that would 
be great.

Thanks,
Jean Yarrington
Independent consultant.
__
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] Counting non-empty levels of a factor

2009-11-08 Thread sylvain willart
Hi everyone,

I'm struggling with a little problem for a while, and I'm wondering if
anyone could help...

I have a dataset (from retailing industry) that indicates which brands
are present in a panel of 500 stores,

store , brand
1 , B1
1 , B2
1 , B3
2 , B1
2 , B3
3 , B2
3 , B3
3 , B4

I would like to know how many brands are present in each store,

I tried:
result <- aggregate(MyData$brand , by=list(MyData$store) , nlevels)

but I got:
Group.1 x
1 , 4
2 , 4
3 , 4

which is not exactly the result I expected
I would like to get sthg like:
Group.1 x
1 , 3
2 , 2
3 , 3

Looking around, I found I can delete empty levels of factor using:
problem.factor <- problem.factor[,drop=TRUE]
But this solution isn't handy for me as I have many stores and should
make a subset of my data for each store before dropping empty factor

I can't either counting the line for each store (N), because the same
brand can appear several times in each store (several products for the
same brand, and/or several weeks of observation)

I used to do this calculation using SAS with:
proc freq data = MyData noprint ; by store ;
 tables  brand / out = result ;
run ;
(the cool thing was I got a database I can merge with MyData)

any idea for doing that in R ?

Thanks in advance,

King Regards,

Sylvain Willart,
PhD Marketing,
IAE Lille, France

__
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] influence.measures(stats): hatvalues(model, ...)

2009-11-08 Thread Peter Ehlers


Sigmund Freud wrote:

Hello:

I am trying to understand the method 'hatvalues(...)', which returns something similar to the diagonals of the plain vanilla hat matrix [X(X'X)^(-1)X'], but not quite. 

A Fortran programmer I am not, but tracing through the code it looks like perhaps some sort of correction based on the notion of 'leave-one-out' variance is being applied. 



I can't see what the problem is. Using the LifeCycleSavings
example from ?influence.measures:

  lm.SR <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)
  X <- model.matrix(lm.SR)
  H <- X %*% solve(t(X) %*% X) %*% t(X)
  hats1 <- diag(H)
  hats2 <- hatvalues(lm.SR)
  all.equal(hats1, hats2)
  #[1] TRUE


Whatever the difference, in simulations 'hatvalues' appears to perform much 
better in the context of identifying outiers using Cook's Distance than the 
diagonals of the plain vanilla hat matrix. (As in 
http://en.wikipedia.org/wiki/Cook's_distance).

Would prefer to understand a little more when using this method. I have 
downloaded the freely available references cited in the help and am in the 
process of digesting them. If someone with knowledge could offer a pointer on 
the most efficient way to get at why 'hatvalues' does what it does, that would 
be great.


In a nutshell, hatvalues are a measure of how unusual
a point is in predictor space, i.e. to what extent it
"sticks out" in one or more of the X-dimensions.

 -Peter Ehlers


Thanks,
Jean Yarrington
Independent consultant.



  
	[[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] lme4 and incomplete block design

2009-11-08 Thread Emmanuel Charpentier
Le dimanche 08 novembre 2009 à 00:05 +0100, Christian Lerch a écrit :
> Dear list members,
> 
> I try to simulate an incomplete block design in which every participants 
> receives 3 out of 4 possible treatment. The outcome in binary.
> 
> Assigning a binary outcome to the BIB or PBIB dataset of the package 
> SASmixed gives the appropriate output.
> With the code below, fixed treatment estimates are not given for each of 
> the 4 possible treatments, instead a kind of summary measure(?) for 
> 'treatment' is given.
> 
> block<-rep(1:24,each=3)
> treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 
> 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 
> 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 
> 4, 4, 1, 3)
> outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
> 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
> 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 
> 0, 0, 1, 0)
> data<-data.frame(block,treatment,outcome)
> lmer(outcome~treatment +(1|block), family=binomial, data=data)
> 
> Is this a problem with the recovery of interblock estimates?

No...

>   Are there 
> special rules how incomplete block designs should look like to enable 
> this recovery?

Neither...

Compare :

> library(lme4)
Le chargement a nécessité le package : Matrix
Le chargement a nécessité le package : lattice
> 
> summary(lmer(outcome~treatment +(1|block), family=binomial,
+  data=data.frame(block<-rep(1:24,each=3),
+treatment<-c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4,
+ 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4,
+ 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4,
+ 2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 1, 3, 2,
+ 4, 4, 1, 3),
+outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 
0,
+   0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 
0,
+   0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 
1,
+   0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 
0))
+  ))
Generalized linear mixed model fit by the Laplace approximation 
Formula: outcome ~ treatment + (1 | block) 
   Data: data.frame(block <- rep(1:24, each = 3), treatment <- c(1, 2,  3, 
2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2,  4, 2, 1, 1, 4, 2, 2, 
3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3, 2,  1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 
3, 3, 1, 2, 4, 2,  1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3), outcome <- c(0, 
0,  0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,  0, 1, 0, 
0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0,  0, 0, 0, 0, 0, 0, 1, 0, 0, 
0, 0, 1, 0, 1, 0, 0, 1, 0, 0,  0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)) 
   AIC  BIC logLik deviance
 86.28 93.1 -40.1480.28
Random effects:
 Groups NameVariance Std.Dev.
 block  (Intercept) 0.60425  0.77734 
Number of obs: 72, groups: block, 24

Fixed effects:
Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.277830.71767  -1.7800.075 .
treatment0.011620.25571   0.0450.964  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
  (Intr)
treatment -0.892

with :

> summary(lmer(outcome~treatment +(1|block), family=binomial,
+  data=data.frame(block<-rep(1:24,each=3),
+treatment<-factor(c(1, 2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 
4,
+2, 4, 3, 3, 1, 2, 4, 2, 1, 1, 4, 2, 2, 3, 
1,
+3, 2, 4, 4, 1, 3, 1, 2, 3, 2, 1, 4, 3, 4, 
1,
+4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1, 2, 4, 2, 
1,
+1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)),
+outcome<-c(0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 
0,
+   0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 
0,
+   0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 
1,
+   0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 
0))
+  ))
Generalized linear mixed model fit by the Laplace approximation 
Formula: outcome ~ treatment + (1 | block) 
   Data: data.frame(block <- rep(1:24, each = 3), treatment <- factor(c(1,  
2, 3, 2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 2, 4, 3, 3, 1,  2, 4, 2, 1, 1, 4, 
2, 2, 3, 1, 3, 2, 4, 4, 1, 3, 1, 2, 3,  2, 1, 4, 3, 4, 1, 4, 3, 2, 1, 3, 4, 
2, 4, 3, 3, 1, 2, 4,  2, 1, 1, 4, 2, 2, 3, 1, 3, 2, 4, 4, 1, 3)), outcome 
<- c(0,  0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,  0, 
0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1,  0, 0, 0, 0, 0, 0, 0, 
1, 0, 0, 0, 0, 1, 0

[R] plot.window arguments being ignored?

2009-11-08 Thread ARRRRR_user

Hi,
Why, when I run the script below, is my y-axis range command being ignored? 
I.e., the y axis graphed consistently fits the line specified in the plot
command, but never fits the line I'm trying to add in the subsequent line
command.
This is my first R script, so if I'm doing anything else wacky that jumps
out, please let me know.  I'm an advanced Stata user, but R is new to me.
The script is designed to simulate a binary-choice decision field theory
example with internally-controlled stop time, as per Neural Networks 19
(2006) 1047-1058.
Thanks in advance!!
-Dan


#-
# Starting preference state; two choices
preference<-c(0,0)
# Preference threshold; when one of the preference states exceeds this
threshhold, deliberation stops
theta<- 2
# Contrast matrix
contrast<-array(c(1, -1, -1, 1), dim=c(2,2))
# Value Matrix; three possible states of the world
value<-array(c(1,0, 0, 1, .5,.5), dim=c(2, 3))
# Feedback Matrix
feedback<-array(c(1, .5, .5, 1), dim=c(2,2))

# Time
t<- array(0, dim=c(1,1))

# Arrays to keep track of history for graphing
time_for_graph<-array(t, dim=c(1,1))
preference_for_graph<-array(preference, dim=c(1,2))

# Moving through time until preference crosses the threshold theta
while (preference[1]http://old.nabble.com/plot.window-arguments-being-ignored--tp26249609p26249609.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] after PCA, the pc values are so large, wrong?

2009-11-08 Thread bbslover

ok,I understand your means, maybe PLS is better for my aim. but I have done
that, also bad. the most questions for me is how to select less variables
from the independent to fit dependent. GA maybe is good way, but I do not
learn it well.

Ben Bolker wrote:
> 
> bbslover  yeah.net> writes:
> 
>> 
> [snip]
> 
>> the fit result below:
>> Call:
>> lm(formula = y ~ x1 + x2 + x3, data = pc)
>> 
>> Residuals:
>>  Min   1Q   Median   3Q  Max 
>> -1.29638 -0.47622  0.01059  0.49268  1.69335 
>> 
>> Coefficients:
>>   Estimate Std. Error t value Pr(>|t|)
>> (Intercept)  5.613e+00  8.143e-02  68.932  < 2e-16 ***
>> x1  -3.089e-05  5.150e-06  -5.998 8.58e-08 ***
>> x2  -4.095e-05  3.448e-05  -1.1880.239
>> x3  -8.106e-05  6.412e-05  -1.2640.210
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
>> 
>> Residual standard error: 0.691 on 68 degrees of freedom
>> Multiple R-squared: 0.3644, Adjusted R-squared: 0.3364 
>> F-statistic: 12.99 on 3 and 68 DF,  p-value: 8.368e-07 
>> 
>> x2,x3 is not significance. by pricipal, after PCA, the pcs should
>> significance, but my data is not, why? 
> 
>   Why is it necessary that the first few principal components
> should have significant relationships with some other response
> values?  The strength, and weakness, of PCA is that it is
> calculated *without regard* to a response variable, so it
> does not constitute "data snooping" ... 
>   I may of course have misinterpreted your question, but at
> a quick look, I don't see anything obviously wrong here.
> 
> __
> 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.
> 
> 

-- 
View this message in context: 
http://old.nabble.com/after-PCA%2C-the-pc-values-are-so-large%2C-wrong--tp26240926p26251658.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] fill map with gradient: package?

2009-11-08 Thread Thomas Steiner
Thanks Paul.
I'm still struggeling with some beginners issues on the ps-import
(windows troubles with installing ghostscript), but when I resolved
them I'm sure that I can use your example code which loos great to me.
Thanks a lot,
Thomas


2009/11/4 Paul Murrell :
> Hi
>
>
> Thomas Steiner wrote:
>>
>> Hi,
>> I'd like to fill an existing svg (or png) map with gradient colors.
>> In detail: The file
>>
>> http://commons.wikimedia.org/wiki/File:Karte_%C3%96sterreich_Bundesl%C3%A4nder.svg
>> should be filled with the population density data from this table:
>> http://de.wikipedia.org/wiki/%C3%96sterreich#Verwaltungsgliederung
>> choosing the right color saturation (or whatever). The final result
>> should be something like this map:
>>
>> http://commons.wikimedia.org/wiki/Image:Bevoelkerungsdichte_-_Oesterreich.svg
>> Is there a package or so for these two tasks (filling and color
>> density ploting)?
>
>
> The 'grImport' package can help with getting the SVG into R (see
> http://www.jstatsoft.org/v30/i04).
>
> First step is to convert the SVG to PostScript (I used InkScape - you can
> play around with how the text comes across, but I'm going to ignore that and
> concentrate on the map regions).
>
> Having done that, the following code loads the map into R and draws it ...
>
> library(grImport)
> PostScriptTrace("Austria-Map-withFonts.ps", charpath=FALSE)
> map <- readPicture("Austria-Map-withFonts.ps.xml")
> grid.picture(map)
>
> ... (the orientation may be 90 degrees out and you may get some warnings
> about character encodings;  the former is easy to fix [see below] and the
> latter can just be ignored for now because we are ignoring the text).  The
> next code shows the breakdown of the map into separate "paths" ...
>
> grid.newpage()
> picturePaths(map)
>
> ... from which we can see that the regions are the first 10 paths ...
>
> grid.newpage()
> grid.picture(map[1:10], use.gc=FALSE)
>
> At this point, you can use grImport to draw the regions with different fill
> colours, or you can just extract the x,y coordinates of the regions and
> go-it-alone.  The following code takes the latter path, setting up 10
> different colours, and drawing each region using grid.polygon().  The
> orientation is fixed by pushing a rotated viewport first ...
>
>
> colours <- hcl(240, 60, seq(30, 80, length=10))
> grid.newpage()
> pushViewport(viewport(angle=-90),
>             grImport:::pictureVP(map[1:10]))
> mapply(function(p, col) {
>           grid.polygon(p$x, p$y, default.units="native",
>                        gp=gpar(fill=col))
>       },
>       regions, colours)
>
>
> Hope that helps.
>
> Paul
>
>
>> Thanks for your help,
>> 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.