date:20090129

Re: [R] Character SNP data to binary MAF data

2009-01-29 Thread Hadassa Brunschwig

Hi

An example is as follows. Consider the character 3x6 matrix:

a A a T A t
G g t T T t
A a C C c c

For each row I would like to identify the most frequent letter and
assign a 1 to it and 0
to the less frequent character. That is, in row 1 the most frequent
letter is A (I do not differentiate between capital and non-capital
letters), in row 2 T and in row 3 C. After the binary conversion
the resulting matrix would look like that:

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

Any suggestions on how to do that (and I am sure I am not the first
one to try this).

Thanks
Hadassa


On Thu, Jan 29, 2009 at 1:50 AM, Jorge Ivan Velez
jorgeivanve...@gmail.com wrote:

 Hi Hadassa,
 Do you have a sample of your data and the output you want? It might be
 useful for us in order to provide any help to you.
 Regards,

 Jorge


 On Wed, Jan 28, 2009 at 8:36 AM, Hadassa Brunschwig
 hadassa.brunsch...@mail.huji.ac.il wrote:

 Hi

 I am sure there is a function out there already but I couldn't find it.
 I have SNP data, that is, a matrix which contains in each row two
 characters (they are different in each row) and I would like to
 convert this matrix to a binary one according to the minor allele
 frequency. For non-geneticists: I want to have a binary matrix
 for which in each row the 0 stands for the less frequent character
 and 1 for the more frequent character.

 Thanks for any suggestions.
 Hadassa

 --
 Hadassa Brunschwig
 PhD Student
 Department of Statistics
 The Hebrew University of Jerusalem
 http://www.stat.huji.ac.il

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





-- 
Hadassa Brunschwig
PhD Student
Department of Statistics
The Hebrew University of Jerusalem
http://www.stat.huji.ac.il

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Odp: Question about collapse/aggregate and avoidance of loops

2009-01-29 Thread Bernd Weiss


markle...@verizon.net schrieb:
 Thanks Petr because I sent Bernd a solution offline but yours is MUCH 
NICER.  it's not worth showing you because it was pretty ugly.




Dear Mark  Petr,

Thank your very much! I like both solutions. Petr's is the more obvious 
one but Mark's solution is good for rethinking my understanding of 
lapply :-)


Again, thank you very much.

Bernd

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Character SNP data to binary MAF data

2009-01-29 Thread Barry Rowlingson

2009/1/29 Hadassa Brunschwig hadassa.brunsch...@mail.huji.ac.il:
 Hi

 An example is as follows. Consider the character 3x6 matrix:

 a A a T A t
 G g t T T t
 A a C C c c

 For each row I would like to identify the most frequent letter and
 assign a 1 to it and 0
 to the less frequent character. That is, in row 1 the most frequent
 letter is A (I do not differentiate between capital and non-capital
 letters), in row 2 T and in row 3 C. After the binary conversion
 the resulting matrix would look like that:

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

 What if there's a tie for most frequent? Do you want 1s for all the
most frequent characters? Or choose one randomly? Or zeroes?

 Examples: what do the following become:

 A A C C T G
 A A C C T T
 A A A A A A

Or are such cases not possible?

 Some hints for you to work on this yourself:

   help('table') - the table function works out counts of elements of vectors
   help('tolower') - for changing upper to lower case
   help('apply') - for working on rows of data frames

 then check out any basic R tutorial on subscripting and replacement,
and you may need to work out how to loop over things with 'for'. You
should be able to make a working solution in a dozen or so lines of R.
Don't be surprised if some R guru on here does it in 2 or 3 lines of
dense, obfuscated stuff!

Barry

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Character SNP data to binary MAF data

2009-01-29 Thread Thomas Lumley



The first step is to convert your data to all uppercase with toupper().

Then it depends on how tidy the data are: are there missing data, are some SNPs 
monomorphic in your sample, etc.

If there are no missing data you can use

N-ncol(the_data)
halfN - N/2

maf_one_row -function(arow) {
   rval-numeric(N)
   if (sum(i-arow==A)halfN) {
rval[]-1
   } else if (sum(i-arow==C)halfN){
rval[i]-1
   } else if (sum(i-arow==T))halfN){
rval[i]-1
   } else if (sum(i-arow==G)halfN){
rval[i]-1
   }
   rval
}

apply(the_data, 1, maf_one_row)

YOu could also use table() to find the two alleles, but you have to make sure 
that the code still works when there is only one allele.

 -thomas

On Thu, 29 Jan 2009, Hadassa Brunschwig wrote:


Hi

An example is as follows. Consider the character 3x6 matrix:

a A a T A t
G g t T T t
A a C C c c

For each row I would like to identify the most frequent letter and
assign a 1 to it and 0
to the less frequent character. That is, in row 1 the most frequent
letter is A (I do not differentiate between capital and non-capital
letters), in row 2 T and in row 3 C. After the binary conversion
the resulting matrix would look like that:

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

Any suggestions on how to do that (and I am sure I am not the first
one to try this).

Thanks
Hadassa


On Thu, Jan 29, 2009 at 1:50 AM, Jorge Ivan Velez
jorgeivanve...@gmail.com wrote:


Hi Hadassa,
Do you have a sample of your data and the output you want? It might be
useful for us in order to provide any help to you.
Regards,

Jorge


On Wed, Jan 28, 2009 at 8:36 AM, Hadassa Brunschwig
hadassa.brunsch...@mail.huji.ac.il wrote:


Hi

I am sure there is a function out there already but I couldn't find it.
I have SNP data, that is, a matrix which contains in each row two
characters (they are different in each row) and I would like to
convert this matrix to a binary one according to the minor allele
frequency. For non-geneticists: I want to have a binary matrix
for which in each row the 0 stands for the less frequent character
and 1 for the more frequent character.

Thanks for any suggestions.
Hadassa

--
Hadassa Brunschwig
PhD Student
Department of Statistics
The Hebrew University of Jerusalem
http://www.stat.huji.ac.il

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







--
Hadassa Brunschwig
PhD Student
Department of Statistics
The Hebrew University of Jerusalem
http://www.stat.huji.ac.il

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



Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Character SNP data to binary MAF data

2009-01-29 Thread Patrick Aboyoun


Hadassa,
You may want to check out the snpMatrix package in Bioconductor

http://bioconductor.org/packages/2.3/bioc/html/snpMatrix.html
http://bioconductor.org/packages/2.4/bioc/html/snpMatrix.html

It contains classes that manage this type of information and should  
minimize your coding effort.



Patrick


Quoting Thomas Lumley tlum...@u.washington.edu:



The first step is to convert your data to all uppercase with toupper().

Then it depends on how tidy the data are: are there missing data, are
some SNPs monomorphic in your sample, etc.

If there are no missing data you can use

N-ncol(the_data)
halfN - N/2

maf_one_row -function(arow) {
   rval-numeric(N)
   if (sum(i-arow==A)halfN) {
rval[]-1
   } else if (sum(i-arow==C)halfN){
rval[i]-1
   } else if (sum(i-arow==T))halfN){
rval[i]-1
   } else if (sum(i-arow==G)halfN){
rval[i]-1
   }
   rval
}

apply(the_data, 1, maf_one_row)

YOu could also use table() to find the two alleles, but you have to
make sure that the code still works when there is only one allele.

 -thomas

On Thu, 29 Jan 2009, Hadassa Brunschwig wrote:


Hi

An example is as follows. Consider the character 3x6 matrix:

a A a T A t
G g t T T t
A a C C c c

For each row I would like to identify the most frequent letter and
assign a 1 to it and 0
to the less frequent character. That is, in row 1 the most frequent
letter is A (I do not differentiate between capital and non-capital
letters), in row 2 T and in row 3 C. After the binary conversion
the resulting matrix would look like that:

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

Any suggestions on how to do that (and I am sure I am not the first
one to try this).

Thanks
Hadassa


On Thu, Jan 29, 2009 at 1:50 AM, Jorge Ivan Velez
jorgeivanve...@gmail.com wrote:


Hi Hadassa,
Do you have a sample of your data and the output you want? It might be
useful for us in order to provide any help to you.
Regards,

Jorge


On Wed, Jan 28, 2009 at 8:36 AM, Hadassa Brunschwig
hadassa.brunsch...@mail.huji.ac.il wrote:


Hi

I am sure there is a function out there already but I couldn't find it.
I have SNP data, that is, a matrix which contains in each row two
characters (they are different in each row) and I would like to
convert this matrix to a binary one according to the minor allele
frequency. For non-geneticists: I want to have a binary matrix
for which in each row the 0 stands for the less frequent character
and 1 for the more frequent character.

Thanks for any suggestions.
Hadassa

--
Hadassa Brunschwig
PhD Student
Department of Statistics
The Hebrew University of Jerusalem
http://www.stat.huji.ac.il

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







--
Hadassa Brunschwig
PhD Student
Department of Statistics
The Hebrew University of Jerusalem
http://www.stat.huji.ac.il

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



Thomas Lumley   Assoc. Professor, Biostatistics
tlum...@u.washington.eduUniversity of Washington, Seattle

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Question about collapse/aggregate and avoidance of loops

2009-01-29 Thread Patrick Burns


I think you are looking for

split(author, id)


Patrick Burns
patr...@burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of The R Inferno and A Guide for the Unwilling S User)

Weiss, Bernd wrote:

Dear all,

given the following data

## original data
id - c(1,1,1,2,2,3)
author - c(A,B,C,D,E,F)
tmp - data.frame(id,author)
tmp


 tmp
  id author
1  1  A
2  1  B
3  1  C
4  2  D
5  2  E
6  3  F

What is the best (most efficient/vectorized/avoiding loops) approach 
to obtain the following data frame?


id author
1A, B, C
2 D, E
3F


Thanks for your help,

Bernd





 version
   _
platform   i386-pc-mingw32
arch   i386
os mingw32
system i386, mingw32
status Patched
major  2
minor  8.1
year   2008
month  12
day22
svn rev47296
language   R
version.string R version 2.8.1 Patched (2008-12-22 r47296)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Character SNP data to binary MAF data

2009-01-29 Thread Barry Rowlingson

2009/1/29 Patrick Aboyoun paboy...@fhcrc.org:
 Hadassa,
 You may want to check out the snpMatrix package in Bioconductor

 http://bioconductor.org/packages/2.3/bioc/html/snpMatrix.html
 http://bioconductor.org/packages/2.4/bioc/html/snpMatrix.html

 It contains classes that manage this type of information and should minimize
 your coding effort.

It's not that much effort - this code turns all ties into 1s:

snp2maf=function(m){
m=toupper(m)
return(t(apply(m,1,makeBin)))
}

makeBin = function(chars){
tc = table(chars)
maxV = names(tc[tc==max(tc)])
matches = match(chars,maxV)
r=as.integer(!is.na(matches))
return(r)
}

 then:

 m
  [,1] [,2] [,3] [,4] [,5]
 [1,] t  g  g  g  t
 [2,] a  G  a  C  c
 [3,] A  T  c  c  C
 [4,] g  T  c  A  C
 [5,] G  C  G  g  G
 [6,] G  t  T  a  C
 [7,] A  G  T  g  T
 [8,] T  a  C  a  T
 [9,] t  g  g  c  T
[10,] A  t  t  c  A
 snp2maf(m)
  [,1] [,2] [,3] [,4] [,5]
 [1,]01110
 [2,]10111
 [3,]00111
 [4,]00101
 [5,]10111
 [6,]01100
 [7,]01111
 [8,]11011
 [9,]11101
[10,]11101


Barry

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Multiple tables

2009-01-29 Thread Gerit Offermann

Dear list,

I have a set of 100+ variables. I would like to have one by one crosstables for 
each variable. I started with
table(variable1, variable2)
table(variable1, variable3)
table(variable1, variable4)
...
table(variable2, variable3)
table(variable2, variable4)
...

It seems rather tedious.

Any better ideas around?

Thanks for any help!
Gerit
-- 
NUR NOCH BIS 31.01.! GMX FreeDSL - Telefonanschluss + DSL 
für nur 16,37 EURO/mtl.!* http://dsl.gmx.de/?ac=OM.AD.PD003K11308T4569a

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Text in a character vector to indicate ifelse argument

2009-01-29 Thread joe1985


Hello

I have a data set that looks like this; 

 b2
  dato chr  status   PRRSvac
PRRSsanVac PRRSsanDk PRRSdk
33  2007-12-03 090432Rød SPF
34  2007-02-09 090432  Rød SPF+sanDK
35  2002-12-17 090432 Rød SPF+DK
36  2002-11-27 090432  Rød SPF+sanDK
37  2002-07-23 090432 Rød SPF+DK
38  2001-08-23 090432Rød SPF
39  2000-01-01 090432  SPF-X,  PRRS-neg.
40  1999-05-01 090432   MS-X,  PRRS-neg.
81  2001-08-23 022458Rød SPF
82  1999-01-22 022458  SPF-X,  PRRS-neg.   
130 2008-10-16 080385 Rød SPF+Myc+Ap2+Nys+DK+Vac   
131 2003-03-18 080385 Rød SPF+Myc+Ap2+DK+Vac
132 2002-11-01 080385 Rød SPF+Myc+DK+Vac
133 2002-02-07 080385Rød SPF+Myc+Vac
134 2000-09-19 080385 MS-X,  PRRS-pos VAC   
135 1999-01-22 080385MS-X,  PRRS-neg
176 2008-10-28 013168 Rød SPF+Myc+Ap2+Nys+DK+Vac
177 2003-05-23 013168 Rød SPF+Myc+Ap2+DK+Vac
178 2002-11-01 013168 Rød SPF+Myc+DK+Vac
179 2001-07-03 013168Rød SPF+Myc+Vac
180 2000-09-01 013168 MS-X,  PRRS-pos VAC  
181 2000-06-02 013168MS-X,  PRRS-neg
182 2000-04-03 013168 SKM-X,  +Ap2,  PRRS-neg   
183 1999-01-22 013168MS-X,  PRRS-neg

Where I have used;

b2$PRRSvac - ifelse(b2$status=='PRRS-pos VAC' | b2$status=='Vac',1,0)
b2$PRRSdk - ifelse(b2$status=='PRRS-pos DK' | b2$status=='DK',1,0)
b2$PRRSsanVac - ifelse(b2$status=='sanVac',1,0)
b2$PRRSsanDk - ifelse(b2$status=='sanDK',1,0)

to creat the last four variables, but it wont work!!! The variable status
has class character. 

Can anyone help me?

-- 
View this message in context: 
http://www.nabble.com/Text-in-a-character-vector-to-indicate-%22ifelse%22-argument-tp21722983p21722983.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] help with plot layout

2009-01-29 Thread Jim Lemon


mau...@alice.it wrote:

It takes a lot of sweat to generate a composite plot with R ...  sigh.
I though I was almost done when I met the umpteenth hurdle. I cannot place a 
nice title on the 2nd plot (raw signal)
on the layout. I do not have control on where either the main option of plot 
function, or title, place the text
string which keeps dysplaying chopped from above. I also tried text, changing many times the string coordinates, but could not see any text anywhere on the canvas . 
By the way, since the layout breaks the canvas into 4 parts, are the text coordinates absolute (referred to the canvas) or 
relative (referred to the part) ?

Please, find attached the generated drawing. The generating script is i the 
following.
  

Hi Maura,
You may find tab.title in the plotrix package helpful.

Jim

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] convergence problem gamm / lme

2009-01-29 Thread geert aarts


Simon, thanks for your reply and your suggestions. 
 
I fitted the following glmm's 
 
gamm3-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=list(code_tripnr=~1),family=poisson))
 
Which worked OK (see summary below)
 
I also fitted a model using quasipoisson, but that didn't help. I actually also 
thought that glmmPQL and gamm estimate the dispersion parameter and hence 
assumes a quasipoisson distribution, even if you specify poisson. Is that 
correct?
 
Finally I tried fitting a model to less data, and sometimes gamm managed to 
converge (see summary below). 
So would it be possible to use the parameter estimates from the model fitted to 
less data as starting values for the gamm fitted to the full data set? 
Or do you have any other suggestions?
 
Thanks.
Cheers Geert
 
 
 
 
 
  
gamm3-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=list(code_tripnr=~1),f
 
amily=poisson))
 
 
 
iteration
1
 
iteration
2
 
iteration
3
 
   detach(Disc_age)
 

summary(gamm3)
 
Linear
mixed-effects model fit by maximum likelihood
 
 Data: NULL
 
  AIC BIC logLik
 
   NA  NA
NA
 
 
 
Random
effects:
 
 Formula: ~1 | code_tripnr
 
(Intercept) Residual
 
StdDev:
0.001391914 231.9744
 
 
 
Variance
function:
 
 Structure: fixed weights
 
 Formula: ~invwt
 
Fixed
effects: count ~ offset(offsetter) + poly(lon, 3) * poly(lat, 3)
 
Value
Std.Error   DF t-value p-value
 
(Intercept)-1.582 11.96 2024 -0.13232174  0.8947
 
poly(lon,
3)1  -4.048   1397.33 2024 -0.00289673  0.9977
 
poly(lon,
3)2 -22.013699.71 2024 -0.03145996  0.9749
 
poly(lon,
3)3  -8.538593.87 2024 -0.01437683  0.9885
 
poly(lat,
3)1-109.624666.05 2024 -0.16458856  0.8693
 
poly(lat,
3)2-104.179381.37 2024 -0.27316977  0.7848
 
poly(lat,
3)3 -10.661221.93 2024 -0.04803585  0.9617
 
poly(lon,
3)1:poly(lat, 3)1  4290.737  61369.98 2024 
0.06991589  0.9443
 
poly(lon,
3)2:poly(lat, 3)1  1853.559  36835.63 2024 
0.05031972  0.9599
 
poly(lon,
3)3:poly(lat, 3)1  -240.521  25771.80 2024 -0.00933272  0.9926
 
poly(lon,
3)1:poly(lat, 3)2  2540.147  41378.38 2024 
0.06138826  0.9511
 
poly(lon,
3)1:poly(lat, 3)2  2540.147  41378.38 2024 
0.06138826  0.9511
 
poly(lon,
3)2:poly(lat, 3)2 -1803.911  21522.17
2024 -0.08381643  0.9332
 
poly(lon,
3)3:poly(lat, 3)2  1040.858  16352.56 2024 
0.06365109  0.9493
 
poly(lon,
3)1:poly(lat, 3)3   632.587  12180.28 2024 
0.05193535  0.9586
 
poly(lon,
3)2:poly(lat, 3)3  -394.339  13088.72 2024 -0.03012818  0.9760
 
poly(lon,
3)3:poly(lat, 3)3  -543.502   6221.71 2024 -0.08735569  0.9304
 
 Correlation:
 
(Intr) ply(ln,3)1
ply(ln,3)2 ply(ln,3)3 ply(lt,3)1
 
poly(lon,
3)10.889
 
poly(lon,
3)20.938  0.878
 
poly(lon,
3)30.843  0.981 
0.792
 
poly(lat,
3)1   -0.829 -0.949 -0.906
-0.882
 
poly(lat,
3)20.859  0.578  0.742 
0.538 -0.474
 
poly(lat,
3)3   -0.552 -0.783 -0.579
-0.756  0.837
 
poly(lon,
3)1:poly(lat, 3)1 -0.947 -0.974
-0.940 -0.940  0.925
 
poly(lon,
3)2:poly(lat, 3)1 -0.934 -0.950
-0.857 -0.929  0.881
 
poly(lon,
3)3:poly(lat, 3)1 -0.818 -0.963
-0.866 -0.945  0.931
 
poly(lon,
3)1:poly(lat, 3)2  0.808  0.975 
0.784  0.968 -0.928
 
poly(lon,
3)2:poly(lat, 3)2  0.737  0.575 
0.853  0.465 -0.659
 
poly(lon,
3)3:poly(lat, 3)2  0.735  0.896 
0.647  0.938 -0.765
 
poly(lon,
3)1:poly(lat, 3)3 -0.794 -0.592
-0.823 -0.518  0.591
 
poly(lon,
3)2:poly(lat, 3)3 -0.542 -0.737
-0.419 -0.781  0.635
 
poly(lon,
3)3:poly(lat, 3)3 -0.398 -0.383
-0.534 -0.334  0.425
 
ply(lt,3)2
ply(lt,3)3 p(,3)1:(,3)1 p(,3)2:(,3)1
 
poly(lon,
3)1
 
poly(lon,
3)2
 
poly(lon,
3)3
 
poly(lat,
3)1
 
poly(lat,
3)2
 
poly(lat,
3)3   -0.136
 
poly(lon,
3)1:poly(lat, 3)1 -0.708  0.690
 
poly(lon,
3)2:poly(lat, 3)1 -0.701  0.710  0.933
 
poly(lon,
3)3:poly(lat, 3)1 -0.499  0.738  0.9560.849
 
poly(lon,
3)1:poly(lat, 3)2  0.458 -0.845
-0.915   -0.934
 
poly(lon,
3)2:poly(lat, 3)2  0.683 -0.344
-0.719   -0.522
 
poly(lon,
3)2:poly(lat, 3)2  0.683 -0.344
-0.719   -0.522
 
poly(lon,
3)3:poly(lat, 3)2  0.464 -0.655
-0.834   -0.884
 
poly(lon,
3)1:poly(lat, 3)3 -0.823  0.241  0.7520.594
 
poly(lon,
3)2:poly(lat, 3)3 -0.300  0.707  0.6120.788
 
poly(lon,
3)3:poly(lat, 3)3 -0.266  0.148  0.4930.250
 
p(,3)3:(,3)1
p(,3)1:(,3)2 p(,3)2:(,3)2 p(,3)3:(,3)2
 
poly(lon,
3)1
 
poly(lon,
3)2
 
poly(lon,
3)3
 
poly(lat,
3)1
 
poly(lat,
3)2
 
poly(lat,
3)3
 
poly(lon,
3)1:poly(lat, 3)1
 
poly(lon,
3)2:poly(lat, 3)1

[R] Odp: Text in a character vector to indicate ifelse argument

2009-01-29 Thread Petr PIKAL

Hi

r-help-boun...@r-project.org napsal dne 29.01.2009 10:25:13:

 
 Hello
 
 I have a data set that looks like this; 
 
  b2
   dato chr  status   PRRSvac
 PRRSsanVac PRRSsanDk PRRSdk
 33  2007-12-03 090432Rød SPF 
 34  2007-02-09 090432  Rød SPF+sanDK 
 35  2002-12-17 090432 Rød SPF+DK 
 36  2002-11-27 090432  Rød SPF+sanDK 
 37  2002-07-23 090432 Rød SPF+DK 
 38  2001-08-23 090432Rød SPF 
 39  2000-01-01 090432  SPF-X,  PRRS-neg. 
 40  1999-05-01 090432   MS-X,  PRRS-neg. 
 81  2001-08-23 022458Rød SPF 
 82  1999-01-22 022458  SPF-X,  PRRS-neg. 
 130 2008-10-16 080385 Rød SPF+Myc+Ap2+Nys+DK+Vac 
 131 2003-03-18 080385 Rød SPF+Myc+Ap2+DK+Vac 
 132 2002-11-01 080385 Rød SPF+Myc+DK+Vac 
 133 2002-02-07 080385Rød SPF+Myc+Vac 
 134 2000-09-19 080385 MS-X,  PRRS-pos VAC 
 135 1999-01-22 080385MS-X,  PRRS-neg 
 176 2008-10-28 013168 Rød SPF+Myc+Ap2+Nys+DK+Vac 
 177 2003-05-23 013168 Rød SPF+Myc+Ap2+DK+Vac 
 178 2002-11-01 013168 Rød SPF+Myc+DK+Vac 
 179 2001-07-03 013168Rød SPF+Myc+Vac 
 180 2000-09-01 013168 MS-X,  PRRS-pos VAC 
 181 2000-06-02 013168MS-X,  PRRS-neg 
 182 2000-04-03 013168 SKM-X,  +Ap2,  PRRS-neg 
 183 1999-01-22 013168MS-X,  PRRS-neg 
 
 Where I have used;
 
 b2$PRRSvac - ifelse(b2$status=='PRRS-pos VAC' | b2$status=='Vac',1,0)
 b2$PRRSdk - ifelse(b2$status=='PRRS-pos DK' | b2$status=='DK',1,0)
 b2$PRRSsanVac - ifelse(b2$status=='sanVac',1,0)
 b2$PRRSsanDk - ifelse(b2$status=='sanDK',1,0)
 
 to creat the last four variables, but it wont work!!! The variable 
status
 has class character. 

What wont work!!! means

 zdrzeni
   sklonot   dobatyp spolfspol.f
1 35  3.00  70.00  stand  35.stand  35.stand
2 50 20.00   9.50  stand  50.stand  50.stand
3 50  5.00  29.50  stand  50.stand  50.stand
4 50 15.00  13.00  stand  50.stand  50.stand


 zdrzeni$v1-ifelse(zdrzeni$typ==stand|zdrzeni$typ==not, 1,0)
 zdrzeni
   sklonot   dobatyp spolfspol.f v1
1 35  3.00  70.00  stand  35.stand  35.stand  1
2 50 20.00   9.50  stand  50.stand  50.stand  1
3 50  5.00  29.50  stand  50.stand  50.stand  1
4 50 15.00  13.00  stand  50.stand  50.stand  1


obviously works, so the problem is in your data and your use of ==.

ifelse(x==abc, 1,0)

means that the result is 1 if x is exactly equivalent to abc and 0 if x 
is  abc, def-abc or in any other variation. You probably need to use 
grep and regexpr expression for test but it is not my cup of tea.

Regards
Petr


 
 Can anyone help me?
 
 -- 
 View this message in context: http://www.nabble.com/Text-in-a-character-
 vector-to-indicate-%22ifelse%22-argument-tp21722983p21722983.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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Graphic device graphics primitives

2009-01-29 Thread Sigbert Klinke


Hi,

I know that some graphics devices in R store graphics primitives such 
that a redraw is possible (e.g. when resizing the window). Is it 
possible to get the current number of stored graphic primitives?


Thanks in advance

Sigbert Klinke

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Graphic device graphics primitives

2009-01-29 Thread Prof Brian Ripley


On Thu, 29 Jan 2009, Sigbert Klinke wrote:


Hi,

I know that some graphics devices in R store graphics primitives such that a 
redraw is possible (e.g. when resizing the window). Is it possible to get the 
current number of stored graphic primitives?


That's not what happens: the graphics engine stores a display list for 
the current device (if enabled), this being a set of graphics calls. 
Devices repaint from either backing store or asking the engine to 
replay the display list.


See ?dev.control. There is no public API to the display list, but of 
course you can interrogate it by C code, with the usual risks.



Thanks in advance

Sigbert Klinke

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


This looks very like a programming question for R-devel: it *is* about 
R's internals.



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

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

[R] Plot dagger symbol in R

2009-01-29 Thread Rau, Roland

Dear all,

I would like to plot the dagger symbol in R (like LaTeX's \dagger).
However, I was unable to do so.

First, I thought maybe dagger actually exists just like the degree
symbol:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(degree))

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(dagger))

However, this was not very successful. New hope emerged that I will
succeed when I read the help page (as so often) for ?plotmath.
There I discovered the 'symbol' thing and read that the Adobe Symbol
font encodings are used. The closest thing I could fine, though, was:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(symbol(\247)))

But this is obviously not a dagger and it seems the Adobe Symbol font
does not have a dagger.

We also know this :-D

library(fortunes)
fortune(Yoda)

So maybe someone can give me some advice?

Thanks in advance,
Roland

--
This mail has been sent through the MPI for Demographic Research.  Should you 
receive a mail that is apparently from a MPI user without this text displayed, 
then the address has most likely been faked. If you are uncertain about the 
validity of this message, please check the mail header or ask your system 
administrator for assistance.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Graphic device graphics primitives

2009-01-29 Thread Duncan Murdoch


Sigbert Klinke wrote:

Hi,

I know that some graphics devices in R store graphics primitives such 
that a redraw is possible (e.g. when resizing the window). Is it 
possible to get the current number of stored graphic primitives?


Thanks in advance

Sigbert Klinke

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
  
You could examine the results of recordPlot, but note the warnings that 
the format is not guaranteed to be stable across R versions.


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.

[R] Use SOCKS proxy

2009-01-29 Thread Daniel Brewer

Hi,

Is there anyway to set up R so that it uses a SOCKS proxy in Linux? I am
getting some strange issues with the Institute's new web filtering
system and want to be able to test whether a problem I have with the
biomaRt package is caused by it.

Many thanks

Dan

-- 
**
Daniel Brewer, Ph.D.

Institute of Cancer Research
Molecular Carcinogenesis
Email: daniel.bre...@icr.ac.uk
**

The Institute of Cancer Research: Royal Cancer Hospital, a charitable Company 
Limited by Guarantee, Registered in England under Company No. 534147 with its 
Registered Office at 123 Old Brompton Road, London SW7 3RP.

This e-mail message is confidential and for use by the a...{{dropped:2}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 dagger symbol in R

2009-01-29 Thread Duncan Murdoch


Rau, Roland wrote:

Dear all,

I would like to plot the dagger symbol in R (like LaTeX's \dagger).
However, I was unable to do so.

First, I thought maybe dagger actually exists just like the degree
symbol:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(degree))

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(dagger))

However, this was not very successful. New hope emerged that I will
succeed when I read the help page (as so often) for ?plotmath.
There I discovered the 'symbol' thing and read that the Adobe Symbol
font encodings are used. The closest thing I could fine, though, was:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(symbol(\247)))

But this is obviously not a dagger and it seems the Adobe Symbol font
does not have a dagger.

We also know this :-D

library(fortunes)
fortune(Yoda)

So maybe someone can give me some advice?
  
I think it will depend on your output device whether the character glyph 
exists, but this works for me on a Mac:


plot(1, main=My dagger: \u2020)

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] Plot dagger symbol in R

2009-01-29 Thread Mark Difford


Hi Roland,

 But this is obviously not a dagger and it seems the Adobe Symbol font
 does not have a dagger.

True, but ... Yoda was here.

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(\u2020))
text(x=0.4, y=0.6, labels=expression(\u2021))

?plotmath, sub Other symbols ... Any Unicode character can be

Regards, Mark.

PS: Works under Windows Vista, but ...


Rau, Roland wrote:
 
 Dear all,
 
 I would like to plot the dagger symbol in R (like LaTeX's \dagger).
 However, I was unable to do so.
 
 First, I thought maybe dagger actually exists just like the degree
 symbol:
 
 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(degree))
 
 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(dagger))
 
 However, this was not very successful. New hope emerged that I will
 succeed when I read the help page (as so often) for ?plotmath.
 There I discovered the 'symbol' thing and read that the Adobe Symbol
 font encodings are used. The closest thing I could fine, though, was:
 
 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(symbol(\247)))
 
 But this is obviously not a dagger and it seems the Adobe Symbol font
 does not have a dagger.
 
 We also know this :-D
 
 library(fortunes)
 fortune(Yoda)
 
 So maybe someone can give me some advice?
 
 Thanks in advance,
 Roland
 
 --
 This mail has been sent through the MPI for Demographic Research.  Should
 you receive a mail that is apparently from a MPI user without this text
 displayed, then the address has most likely been faked. If you are
 uncertain about the validity of this message, please check the mail header
 or ask your system administrator for assistance.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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://www.nabble.com/Plot-dagger-symbol-in-R-tp21725141p21725439.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] Plot dagger symbol in R

2009-01-29 Thread baptiste auguie


Hi,

If all else fails, you could consider using LaTeX itself with psfrag,  
or perhaps a similar idea involving eps2pgf.


http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/PsFrag

Hope this helps,

baptiste

On 29 Jan 2009, at 11:24, Rau, Roland wrote:


Dear all,

I would like to plot the dagger symbol in R (like LaTeX's \dagger).
However, I was unable to do so.

First, I thought maybe dagger actually exists just like the degree
symbol:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(degree))

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(dagger))

However, this was not very successful. New hope emerged that I will
succeed when I read the help page (as so often) for ?plotmath.
There I discovered the 'symbol' thing and read that the Adobe Symbol
font encodings are used. The closest thing I could fine, though, was:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(symbol(\247)))

But this is obviously not a dagger and it seems the Adobe Symbol font
does not have a dagger.

We also know this :-D

library(fortunes)
fortune(Yoda)

So maybe someone can give me some advice?

Thanks in advance,
Roland

--
This mail has been sent through the MPI for Demographic Research.   
Should you receive a mail that is apparently from a MPI user without  
this text displayed, then the address has most likely been faked. If  
you are uncertain about the validity of this message, please check  
the mail header or ask your system administrator for assistance.


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


_

Baptiste Auguié

School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 dagger symbol in R

2009-01-29 Thread Prof Brian Ripley


On Thu, 29 Jan 2009, Mark Difford wrote:



Hi Roland,


But this is obviously not a dagger and it seems the Adobe Symbol font
does not have a dagger.


True, but ... Yoda was here.

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(\u2020))
text(x=0.4, y=0.6, labels=expression(\u2021))

?plotmath, sub Other symbols ... Any Unicode character can be

Regards, Mark.

PS: Works under Windows Vista, but ...


That would be expected to work

(in a UTF-8 locale || on Windows)
 on a device with Unicode support
 in a font that has the glyph.

(it works on Windows because we fake much of a UTF-8 locale there).

There is an alternative: dagger _is_ in the standard Adobe character 
set, so this will work as something \206 (untested) in 8-bit Windows 
character sets on windows(), postscript(), pdf() ... devices


As ever, the information asked for in the posting guide helps us give 
a better answer.





Rau, Roland wrote:


Dear all,

I would like to plot the dagger symbol in R (like LaTeX's \dagger).
However, I was unable to do so.

First, I thought maybe dagger actually exists just like the degree
symbol:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(degree))

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(dagger))

However, this was not very successful. New hope emerged that I will
succeed when I read the help page (as so often) for ?plotmath.
There I discovered the 'symbol' thing and read that the Adobe Symbol
font encodings are used. The closest thing I could fine, though, was:

plot(0:1,0:1, type=n)
text(x=0.5, y=0.5, labels=expression(symbol(\247)))

But this is obviously not a dagger and it seems the Adobe Symbol font
does not have a dagger.


But the Adobe Standard encoding does.


We also know this :-D

library(fortunes)
fortune(Yoda)

So maybe someone can give me some advice?

Thanks in advance,
Roland


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

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

Re: [R] Multiple tables

2009-01-29 Thread jim holtman

you can use combn to create the combinations and the following will
create a list of all the results:

x1 - x2 - x3 - x4 - 1:10
comb - combn(c('x1','x2', 'x3', 'x4'), 2)
myTab - lapply(seq(ncol(comb)), function(x){
table(get(comb[1, x]), get(comb[2, x]))
})
# put names of the combinations
names(myTab) - apply(comb, 2, paste, collapse=:)



On Thu, Jan 29, 2009 at 4:19 AM, Gerit Offermann gerit.offerm...@gmx.de wrote:
 Dear list,

 I have a set of 100+ variables. I would like to have one by one crosstables 
 for each variable. I started with
 table(variable1, variable2)
 table(variable1, variable3)
 table(variable1, variable4)
 ...
 table(variable2, variable3)
 table(variable2, variable4)
 ...

 It seems rather tedious.

 Any better ideas around?

 Thanks for any help!
 Gerit
 --
 NUR NOCH BIS 31.01.! GMX FreeDSL - Telefonanschluss + DSL
 für nur 16,37 EURO/mtl.!* http://dsl.gmx.de/?ac=OM.AD.PD003K11308T4569a

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] ANOVA in R

2009-01-29 Thread Amit Patel

Hi

I Have a very large dataset that I would like to conduct ANOVA tests on. Im not 
a very strong programmer so any help would be appreciated.
the format is 

IdentifierÂ Â Â Â Â Â Â Â Â Â Â Â  A1Â Â Â Â Â Â  A2Â Â Â Â Â Â Â  B1Â Â Â Â Â  
B2Â Â Â Â Â Â  C1Â Â  C2Â Â Â Â Â  Norm1Â Â Â Â Â Â Â Â  Norm2
1234Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â  1Â Â Â Â Â Â Â  1Â Â Â Â Â Â Â Â Â Â Â  
NAÂ Â Â Â  NAÂ Â Â Â Â  4Â Â Â Â Â Â  3Â Â Â Â Â Â Â  
NAÂ Â Â Â Â Â Â Â Â Â Â Â Â Â  NA
4567Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â  2Â Â Â Â Â Â Â  
2Â Â Â Â Â Â Â Â Â Â Â Â Â  4Â Â Â Â Â  4Â Â Â Â Â Â Â Â  8Â Â Â Â Â Â  
8Â Â Â Â Â Â  9Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â  9


and so on
I have 10 runs for 3 different doses plus the normal state. Any help greatly 
appreciated

 
 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 
 
  
  

  

  

  

  

  

 





  
[[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] Plot dagger symbol in R

2009-01-29 Thread Rau, Roland

Dear all,

thank you very much.
Yes, indeed, I should have included the sessionInfo().
Both examples (unicode and standard encoding) work fine; however, I
prefer the suggestion by Prof. Ripley since it allows me to do things
like:

plot(1, main=expression(e^\206))

This fails, however, with:

plot(1, main=expression(e^\u2020))

Thank you once again,
Roland

 sessionInfo()
R version 2.7.0 (2008-04-22) 
i386-pc-mingw32 

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

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

 


 -Original Message-
 From: r-help-boun...@r-project.org 
 [mailto:r-help-boun...@r-project.org] On Behalf Of Prof Brian Ripley
 Sent: Thursday, January 29, 2009 1:35 PM
 To: Mark Difford
 Cc: r-help@r-project.org
 Subject: Re: [R] Plot dagger symbol in R
 
 On Thu, 29 Jan 2009, Mark Difford wrote:
 
 
  Hi Roland,
 
  But this is obviously not a dagger and it seems the Adobe 
 Symbol font
  does not have a dagger.
 
  True, but ... Yoda was here.
 
  plot(0:1,0:1, type=n)
  text(x=0.5, y=0.5, labels=expression(\u2020))
  text(x=0.4, y=0.6, labels=expression(\u2021))
 
  ?plotmath, sub Other symbols ... Any Unicode character 
 can be
 
  Regards, Mark.
 
  PS: Works under Windows Vista, but ...
 
 That would be expected to work
 
 (in a UTF-8 locale || on Windows)
  on a device with Unicode support
  in a font that has the glyph.
 
 (it works on Windows because we fake much of a UTF-8 locale there).
 
 There is an alternative: dagger _is_ in the standard Adobe character 
 set, so this will work as something \206 (untested) in 8-bit Windows 
 character sets on windows(), postscript(), pdf() ... devices
 
 As ever, the information asked for in the posting guide helps us give 
 a better answer.
 
 
 
  Rau, Roland wrote:
 
  Dear all,
 
  I would like to plot the dagger symbol in R (like LaTeX's \dagger).
  However, I was unable to do so.
 
  First, I thought maybe dagger actually exists just like the degree
  symbol:
 
  plot(0:1,0:1, type=n)
  text(x=0.5, y=0.5, labels=expression(degree))
 
  plot(0:1,0:1, type=n)
  text(x=0.5, y=0.5, labels=expression(dagger))
 
  However, this was not very successful. New hope emerged that I will
  succeed when I read the help page (as so often) for ?plotmath.
  There I discovered the 'symbol' thing and read that the 
 Adobe Symbol
  font encodings are used. The closest thing I could fine, 
 though, was:
 
  plot(0:1,0:1, type=n)
  text(x=0.5, y=0.5, labels=expression(symbol(\247)))
 
  But this is obviously not a dagger and it seems the Adobe 
 Symbol font
  does not have a dagger.
 
 But the Adobe Standard encoding does.
 
  We also know this :-D
 
  library(fortunes)
  fortune(Yoda)
 
  So maybe someone can give me some advice?
 
  Thanks in advance,
  Roland
 
 -- 
 Brian D. Ripley,  rip...@stats.ox.ac.uk
 Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
 University of Oxford, Tel:  +44 1865 272861 (self)
 1 South Parks Road, +44 1865 272866 (PA)
 Oxford OX1 3TG, UKFax:  +44 1865 272595
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 

--
This mail has been sent through the MPI for Demographic Research.  Should you 
receive a mail that is apparently from a MPI user without this text displayed, 
then the address has most likely been faked. If you are uncertain about the 
validity of this message, please check the mail header or ask your system 
administrator for assistance.

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

[R] t-test

2009-01-29 Thread Amit Patel

When doing the t-test in the below manner will r compare each element of the 
array with the relevant one. I.e. if i was comparing x and y would (1 and 0) 
and (1 and 9) be treated as separate variables. Or does it just assume one 
variable. 


# test data

x - c(1,1.1,1.15,1.2,1.21,1.23)

y - c(0.9,1,1.16,1.18,1.19,1.2)

z - c(1.4,1.42,1.43,1.44,1.45,1.46)


###Â  Student's t-test

# for help in R type ?t.test()

# defaults are:
# alternative = two.sided i.e. two-sided t-test
# var.equal = FALSE i.e. unequal variance

# note:
# na.rm = TRUE removes missing values
# $p.value gives the p-value for the test

t.test(x, y, na.rm=TRUE, paired=FALSE)$p.value

# gives 0.5026467

t.test(x, z, na.rm=TRUE, paired=FALSE)$p.value

# gives 0.0003166352



  
[[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] Goodness of fit for gamma distributions

2009-01-29 Thread Dan31415


Ah yes, that does produce a nice plot. Can i just ask what exactly it is
showing. It seems to me to be a sort of Q-Q plot but with a different set of
axes. Is this correct, if so do the same interpretation rules apply for this
plot, i.e. departures from either end of the curve show poor fitting of the
extreme data.

thanks for your help Remko, its been very helpful.

Dann



Remko Duursma-2 wrote:
 
 It sounds like you just want to graph it though. For gammas, it's nice
 to graph the log of the density, because
 the tail is so thin and long, so you don't see much otherwise:
 
 mydata - rgamma(1, shape=1.1, rate=2.5)
 
 # now suppose you fit a gamma distribution, and get these estimated
 parameters:
 shapeest - 1.101
 rateest - 2.49
 
 h - hist(mydata, breaks=50, plot=FALSE)
 plot(h$mids, log(h$density))
 curve(log(dgamma(x, shape=shapeest, rate=rateest)), add=TRUE)
 
 
 #Remko
 
 
 -
 Remko Duursma
 Post-Doctoral Fellow
 
 Centre for Plant and Food Science
 University of Western Sydney
 Hawkesbury Campus
 Richmond NSW 2753
 
 Dept of Biological Science
 Macquarie University
 North Ryde NSW 2109
 Australia
 
 Mobile: +61 (0)422 096908
 
 
 
 On Wed, Jan 28, 2009 at 1:13 AM, Dan31415 d.m.mitch...@reading.ac.uk
 wrote:

 Thanks for that Remko, but im slightly confused because isnt this testing
 the
 goodness of fit of 2 slightly different gamma distributions, not of how
 well
 a gamma distribution is representing the data.

 e.g.

 data.vec-as.vector(data)

 (do some mle to find the parameters of a gamma distribution for data.vec)

 xrarea-seq(-2,9,0.05)
 yrarea-dgamma(xrarea,shape=7.9862,rate=2.6621)

 so now yrarea is the gamma distribution and i want to compare it with
 data.vec to see how well it fits.

 regards,
 Dann


 Remko Duursma-2 wrote:

 Hi Dann,

 there is probably a better way to do this, but this works anyway:

 # your data
 gamdat - rgamma(1, shape=1, rate=0.5)

 # comparison to gamma:
 gamsam - rgamma(1, shape=1, rate=0.6)

 qqplot(gamsam,gamdat)
 abline(0,1)


 greetings
 Remko


 -
 Remko Duursma
 Post-Doctoral Fellow

 Centre for Plant and Food Science
 University of Western Sydney
 Hawkesbury Campus
 Richmond NSW 2753

 Dept of Biological Science
 Macquarie University
 North Ryde NSW 2109
 Australia

 Mobile: +61 (0)422 096908



 On Tue, Jan 27, 2009 at 3:38 AM, Dan31415 d.m.mitch...@reading.ac.uk
 wrote:

 I'm looking for goodness of fit tests for gamma distributions with
 large
 data
 sizes. I have a matrix with around 10,000 data values in it and i have
 fitted a gamma distribution over a histogram of the data.

 The problem is testing how well that distribution fits. Chi-squared
 seems
 to
 be used more for discrete distributions and kolmogorov-smirnov seems
 that
 large sample sizes make it had to evaluate the D statistic. Also i
 haven't
 found a qq plot for gamma, although i think this might be an
 appropriate
 test.

 in summary
 -is there a gamma goodness of fit test that doesnt depend on the sample
 size?
 -is there a way of using qqplot for gamma distributions, if so how
 would
 you
 calculate it from a matrix of data values?

 regards,
 Dann
 --
 View this message in context:
 http://www.nabble.com/Goodness-of-fit-for-gamma-distributions-tp21668711p21668711.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-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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://www.nabble.com/Goodness-of-fit-for-gamma-distributions-tp21668711p21686095.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-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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://www.nabble.com/Goodness-of-fit-for-gamma-distributions-tp21668711p21725468.html
Sent from the R help mailing list archive at Nabble.com.

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

[R] Inconsistency in F values from dropterm and anova

2009-01-29 Thread IMS


Hi,

I'm working on fitting a glm model to my data using Gamma error structure
and reciprocal link.  I've been using dropterm (MASS) in the model
simplification process, but the F values from analysis of deviance tables
reported by dropterm and anova functions are different - sometimes
significantly so.  However, the reported residual deviances, degrees of
freedom, etc. are not different.  I don't know how to calculate F values
from deviance tables (and can't seem to find anything suggesting how), and
so haven't been able to calculate F for myself to see which is more
accurate.  Below is an example of the code, and I'm using R version 2.8.0.

 model1=glm(diff2^(0.491)~mtype*morder,family=Gamma)
 dropterm(model1,test=F)
Single term deletions

Model:
diff2^(0.491) ~ mtype * morder
Df  Deviance AIC   F value   Pr(F)
   23.181   -16.813   
mtype:morder  3   24.729   -13.741  2.0694   0.1096

 model2=update(model1,~.-mtype:morder)
 anova(model1,model2,test=F)
Analysis of Deviance Table

Model 1: diff2^(0.491) ~ mtype * morder
Model 2: diff2^(0.491) ~ mtype + morder
  Resid. Df Resid. Dev Df Deviance  F  Pr(F)  
19323.1814 
29624.7288 -3  -1.5475 3.0241 0.03352 *
---
Signif. codes:  0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â 
â 1 

-- 
View this message in context: 
http://www.nabble.com/Inconsistency-in-F-values-from-dropterm-and-anova-tp21725668p21725668.html
Sent from the R help mailing list archive at Nabble.com.

[[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] Dynamic random effects model

2009-01-29 Thread Ben Bolker

Joseph Magagnoli jcm331 at gmail.com writes:

 
 All R experts,
 How do I fit a dynamic Random effects model with a binary dependent variable
 in R
 Thanks
 JCM
 

  You haven't given us nearly enough information to go on.
If you're talking about something like a state-space model with
a binary response, I would probably say your best bet is
a Bayesian approach, prob. via JAGS/R2jags or WinBUGS/R2WinBUGS.
(A lot) more context will give a higher probability of a useful
answer.

  good luck
Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 dagger symbol in R

2009-01-29 Thread Duncan Murdoch


On 1/29/2009 8:56 AM, Rau, Roland wrote:

Dear all,

thank you very much.
Yes, indeed, I should have included the sessionInfo().
Both examples (unicode and standard encoding) work fine; however, I
prefer the suggestion by Prof. Ripley since it allows me to do things
like:

plot(1, main=expression(e^\206))

This fails, however, with:

plot(1, main=expression(e^\u2020))  

Thank you once again,
Roland


sessionInfo()
R version 2.7.0 (2008-04-22) 


That's 9 months old  -- ancient!  The exact line you posted works for me 
in 2.8.1.  However, I don't like the dagger in the default font; it 
looks better as


plot(1, main=expression(e^\u2020), family=serif)


Duncan Murdoch

i386-pc-mingw32 


locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

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







-Original Message-
From: r-help-boun...@r-project.org 
[mailto:r-help-boun...@r-project.org] On Behalf Of Prof Brian Ripley

Sent: Thursday, January 29, 2009 1:35 PM
To: Mark Difford
Cc: r-help@r-project.org
Subject: Re: [R] Plot dagger symbol in R

On Thu, 29 Jan 2009, Mark Difford wrote:


 Hi Roland,

 But this is obviously not a dagger and it seems the Adobe 
Symbol font

 does not have a dagger.

 True, but ... Yoda was here.

 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(\u2020))
 text(x=0.4, y=0.6, labels=expression(\u2021))

 ?plotmath, sub Other symbols ... Any Unicode character 
can be


 Regards, Mark.

 PS: Works under Windows Vista, but ...

That would be expected to work

(in a UTF-8 locale || on Windows)
 on a device with Unicode support
 in a font that has the glyph.

(it works on Windows because we fake much of a UTF-8 locale there).

There is an alternative: dagger _is_ in the standard Adobe character 
set, so this will work as something \206 (untested) in 8-bit Windows 
character sets on windows(), postscript(), pdf() ... devices


As ever, the information asked for in the posting guide helps us give 
a better answer.




 Rau, Roland wrote:

 Dear all,

 I would like to plot the dagger symbol in R (like LaTeX's \dagger).
 However, I was unable to do so.

 First, I thought maybe dagger actually exists just like the degree
 symbol:

 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(degree))

 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(dagger))

 However, this was not very successful. New hope emerged that I will
 succeed when I read the help page (as so often) for ?plotmath.
 There I discovered the 'symbol' thing and read that the 
Adobe Symbol
 font encodings are used. The closest thing I could fine, 
though, was:


 plot(0:1,0:1, type=n)
 text(x=0.5, y=0.5, labels=expression(symbol(\247)))

 But this is obviously not a dagger and it seems the Adobe 
Symbol font

 does not have a dagger.

But the Adobe Standard encoding does.

 We also know this :-D

 library(fortunes)
 fortune(Yoda)

 So maybe someone can give me some advice?

 Thanks in advance,
 Roland

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

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

and provide commented, minimal, self-contained, reproducible code.



--
This mail has been sent through the MPI for Demographic Research.  Should you 
receive a mail that is apparently from a MPI user without this text displayed, 
then the address has most likely been faked. If you are uncertain about the 
validity of this message, please check the mail header or ask your system 
administrator for assistance.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Question On CrossTable function in gmodels package

2009-01-29 Thread eugen pircalabelu

Hi R-users,
I have the following problem with CrossTable function within “gmodels” package: 
the output of the function (format “spss” and asresid=T) can not be stored 
within another object.
For example:

library(gmodels)
data(infert, package = datasets)
 CrossTable(infert$education, infert$induced)-aa # the function prints 
 everything ok on the screen
 aa  # works just fine
$t
 y
x  0  1  2
  0-5yrs   4  2  6
  6-11yrs 78 27 15
  12+ yrs 61 39 16
….
But when I call
CrossTable(infert$education, infert$induced, asresid=TRUE, format=SPSS)-aaa 
 # the function prints everything ok on the screen
 aaa  # now I have a problem 
NULL

Why is aaa object NULL? Should it be NULL? I suspect it has something to do 
with the SPSS-format option but I’m not sure 
The reason I want it stored in an object, is to use the adjusted standardized 
residuals in a LateX document with xtable. 
Does anyone have a solution, please?

Thank you very much and have a great day ahead!  


PS: 
sessionInfo()
R version 2.7.1 (2008-06-23) 
i386-pc-mingw32 

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United 
States.1252;LC_MONETARY=English_United 
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

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

other attached packages:
[1] gmodels_2.14.1 xtable_1.5-4  

loaded via a namespace (and not attached):
[1] gdata_2.3.1  gtools_2.4.0 MASS_7.2-42  tools_2.7.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.

Re: [R] t-test

2009-01-29 Thread patricia garcía gonzález


Hi, 

If you mean if the t.test is done as if the samples where paired, the answer is 
yes if you write argument paired = TRUE; and the pairs are done in order, that 
is 1º with 1º, 2º with 2º, etc.

As you wrote, (paired = FALSE) the t.test is unpaired, and the order of 
elements in the vectors are not taken into account.

Regards

Patricia


 Date: Thu, 29 Jan 2009 13:59:41 +
 From: amitrh...@yahoo.co.uk
 To: r-help@r-project.org
 Subject: [R] t-test
 
 When doing the t-test in the below manner will r compare each element of the 
 array with the relevant one. I.e. if i was comparing x and y would (1 and 0) 
 and (1 and 9) be treated as separate variables. Or does it just assume one 
 variable. 
 
 
 # test data
 
 x - c(1,1.1,1.15,1.2,1.21,1.23)
 
 y - c(0.9,1,1.16,1.18,1.19,1.2)
 
 z - c(1.4,1.42,1.43,1.44,1.45,1.46)
 
 
 ###Â  Student's t-test
 
 # for help in R type ?t.test()
 
 # defaults are:
 # alternative = two.sided i.e. two-sided t-test
 # var.equal = FALSE i.e. unequal variance
 
 # note:
 # na.rm = TRUE removes missing values
 # $p.value gives the p-value for the test
 
 t.test(x, y, na.rm=TRUE, paired=FALSE)$p.value
 
 # gives 0.5026467
 
 t.test(x, z, na.rm=TRUE, paired=FALSE)$p.value
 
 # gives 0.0003166352
 
 
 
   
   [[alternative HTML version deleted]]
 

_


[[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] Repeated measures design for GAM? - corrected question...

2009-01-29 Thread Simon Wood

Diederik,

The `gamm' call looks fine provided `park' is nested in `count'. You are 
basically adding a random effect for each survey, and a random effect for 
each `park' within survey to the the smooth effects. 

best,
Simon

ps. I assume there are more than 18 observations in the full dataset --- the 
model sturcture is a bit too complicated otherwise...

On Wednesday 28 January 2009 19:44, Strubbe Diederik wrote:
 Dear Simon,
 Many thanks for pointing me to the GAMM! For clarification, Bird_abundance
 are breeding densities ( e.g. 1.25 BP/ha, 2.20 BP/ha,...) and count is just
 the actual survey(e.g. first_survey,...). The dataset looks like
 Bird_abundanceStudy_area  YEARCOUNT   X1  X2  X3  
 X4
 0.15  area_1  2004first_survey�   �   �   �
 1.26  area_1  2004second_survey   �   �   �   �
 2.47  area_1  2005third_survey�   �   �   �
 0.00  area_1  2005fourth_survey   �   �   �   �
 0.23  area_1  2006fifht_survey�   �   �   �
 2.64  area_1  2006sixth_survey�   �   �   �
 4.14  area_2  2004first_survey�   �   �   �
 5.00  area_2  2004second_survey   �   �   �   �
 6.80  area_2  2005third_survey�   �   �   �
 0.15  area_2  2005fourth_survey   �   �   �   �
 0.25  area_2  2006fifht_survey�   �   �   �
 2.36  area_2  2006sixth_survey�   �   �   �
 2.59  area_3  2004first_survey�   �   �   �
 6.31  area_3  2004second_survey   �   �   �   �
 0.15  area_3  2005third_survey�   �   �   �
 2.85  area_3  2005fourth_survey   �   �   �   �
 2.48  area_3  2006fifht_survey�   �   �   �
 1.23  area_3  2006sixth_survey�   �   �   �
 � �   �   �   �   �   �   �

 Am I correct in assuming the following is a valid syntax for this repeated
 measures design?:

 model -gamm(Bird_abundance ~ YEAR + s(X1)+ s(X2)+ s(X3)+
 s(X4),random=list(count=~1,park=~1))

 best wishes and thanks again,

 Diederik





 Diederik Strubbe
 Evolutionary Ecology Group
 Department of Biology, University of Antwerp
 Universiteitsplein 1
 B-2610 Antwerp, Belgium
 http://webhost.ua.ac.be/deco
 tel : 32 3 820 23 85



 -Original Message-
 From: Strubbe Diederik
 Sent: Wed 28-1-2009 18:09
 To: r-help@R-project.org
 Subject: Repeated measures design for GAM? - corrected question...

 Dear all,

 I have a question on the use of GAM with repeated measures. My dataset is
 as follows: - a number of study areas where bird abundance has been
 determined. Counts have been performed in 3 consecutive years and there
 were 2 counts per year (i.e. in total 6 counts). - a number of
 environmental predictors that do not change over year Xi). When using a
 GLM, a repeated measures design would like: (for example)

 lme(Bird_abundance = study_area + count +year+ X1 + X2 + X3,random =
 ~count|study_area).

 However, I have found no analogue design for a GAM. For now, I have
 averaged my bird abundances but I wondered whether a more subtle and
 elegant strategy exists...?

 Many thanks,


 Diederik

 Diederik Strubbe
 Evolutionary Ecology Group
 Department of Biology, University of Antwerp
 Universiteitsplein 1
 B-2610 Antwerp, Belgium
 http://webhost.ua.ac.be/deco
 tel : 32 3 820 23 85




   [[alternative HTML version deleted]]

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Stacking data

2009-01-29 Thread Amit Patel

Hi

I have data in the format below

Â Â Â Â  AgeÂ Â  Â  Â  Â  Â Â  V1Â Â Â Â Â Â  V2 Â  Â Â  V3Â Â Â Â Â Â  V4
Â Â  23646Â Â Â Â Â Â Â Â  45190 50333 55166 56271
Â Â  26174Â Â Â Â Â Â Â Â  35535 38227 37911 41184
Â Â  27723Â Â Â Â Â Â Â Â Â  25691 25712 26144 26398

and would like to sort it as follows

Â Â Â Â Â Â Â Â Â Â Â Â  Â  Â  AgeÂ Â Â Â Â  values Â  Â  ind
Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â  23646Â Â Â  45190 Â Â  V1
Â Â Â Â Â Â Â Â Â Â Â Â  Â  26174Â Â  35535 Â  Â Â  V1
Â Â Â Â Â Â Â Â Â Â Â Â Â Â  27723Â Â Â  25691 Â  Â  V2
Â Â Â Â Â Â Â Â Â Â  Â Â  27193Â Â Â  30949 Â  Â Â  V2

But i have 41 columns (age columnÂ  +Â  40 individuals) 
I have the following script but an error is thrown up
can anyone help, where am i going wrong

zz - read.csv(Filename.csv,strip.white = TRUE)
zzz - cbind(zz[gl(nrow(zz), 1, 40*nrow(zz)), 1], stack(zz[, 2:41]))

ERROR MESSAGE
Error in data.frame(..., check.names = FALSE) : 
Â  arguments imply differing number of rows: 3789080, 0




  
[[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] ggplot2 - how to change location / position of wind rose axis labels?

2009-01-29 Thread hadley wickham

Hi Darren,

Sorry for taking so long to get back to you - my hard drive died and
it's taking me a while to get up and running again of backups etc.

One solution to your problem is to draw the axis labels yourself:

labels - data.frame(rrating = 1:10, count = 19000)

awind +
   opts(axis.text.x = theme_blank()) +
   geom_text(aes(label = rrating, y = count, fill = NULL), data = labels)

Adjusting count as needed to get them in the right position.  I'll
think about how to do this better in general.

Regards,

Hadley

On Fri, Jan 23, 2009 at 4:39 PM, Darren Norris doo...@hotmail.com wrote:

 Dear R users,
 First just want to say thank you to all for developing such a wonderful
 software and packages.

 I need to produce a wind rose plot. Tried with packages circular and plotrix
 and couldn't quite get what I want. Moved to package ggplot2 and it's going
 great. However stuck in how to move axis labels.

 I am using the wind rose from the help to learn how to do what I need (code
 example below). The plot produced has axis labels on top of the axis line
 (not sure if it's actually the plot or axis line...) - which means the
 labels are unclear for what I need to produce.

 How can I move the labels to be outside of the line? I have read the online
 book ( http://had.co.nz/ggplot2/book/ ), help, tried changing various theme
 and scale settings and searched the package website and mailing list (
 http://had.co.nz/ggplot2/ ). If the answer is there I'm too stupid to see it
 - do I need to play with grobs? If so how?

 Any help much appreciated,
 Darren

 R version 2.8.1 (2008-12-22)
 i386-pc-mingw32

 locale:
 LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
 Kingdom.1252;LC_MONETARY=English_United
 Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

 other attached packages:
 [1] ggplot2_0.8.1 reshape_0.8.2 plyr_0.1.4proto_0.3-8


 #Using an example form the coord_polar help
 library(ggplot2)
 movies$rrating - factor(round_any(movies$rating, 1))
 movies$budgetq - factor(chop(movies$budget, 4), labels = 1:4)
 doh - ggplot(movies, aes(x = rrating, fill = budgetq))
 doh + geom_bar(width = 1) + coord_polar()

 #Now with my theme (hacked from theme_bw) getting close to what I need

 theme_bwdn-function (base_size = 12)
 {
structure(list(axis.line = theme_blank(), axis.text.x = theme_text(size
 = base_size *
   1, lineheight = 0.9, vjust = 1), axis.text.y = theme_blank(),
 axis.ticks = theme_segment(colour = black,
size = 0.2), axis.title.x = theme_blank(), axis.title.y =
 theme_blank(), axis.ticks.length = unit(0,
lines), axis.ticks.margin = unit(0, lines), legend.background =
 theme_rect(colour = NA),
legend.key = theme_rect(colour = grey80), legend.key.size =
 unit(1.2,
lines), legend.text = theme_text(size = base_size *
0.8), legend.title = theme_text(size = base_size *
0.8, face = bold, hjust = 0), legend.position = right,
panel.background = theme_rect(fill = white, colour = NA),
panel.border = theme_rect(fill = NA, colour = grey50),
panel.grid.major = theme_line(colour = black, size = 0.2),
panel.grid.minor = theme_line(colour = black, size = 0.5),
panel.margin = unit(0.25, lines), strip.background =
 theme_rect(fill = grey80,
colour = grey50), strip.label = function(variable,
value) value, strip.text.x = theme_text(size = base_size *
0.8), strip.text.y = theme_text(size = base_size *
0.8, angle = -90), plot.background = theme_rect(colour = NA),
plot.title = theme_text(size = base_size * 1.2), plot.margin =
 unit(rep(0.5,
4), lines)), class = options)
 }

 awind-doh + geom_bar(width = 1) + coord_polar()

 #produces the wind rose but how to move axis labels?
 awind + theme_bwdn()




 --
 View this message in context: 
 http://www.nabble.com/ggplot2---how-to-change-location---position-of-wind-rose-axis-labels--tp21635526p21635526.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.




-- 
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] convergence problem gamm / lme

2009-01-29 Thread Simon Wood

Geert,

Sorry for slow reply... I don't see any obvious problems with what you've 
done, so I guess it's the usual problem that PQL just doesn't *have* to 
converge, and the bit of extra flexibility of using a smooth is too much for 
it in this case. If you send me the data offline I can dig a little bit more 
if you like (I'll only use the data for this purpose etc. etc.) 

You are right that PQL does the same thing for Poisson and quasi-poisson. I 
don't think there is an easy way to use the values for the reduced dataset 
fit in the full dataset fitting, unfortunately. 

Another option is to use `gam' to fit the random effects. It'll be a bit slow 
with 70+ random effects, as you have, and it's a bit more work to set up, but 
it should converge. See ?gam.models which has some examples showing how to do 
this.

best,
Simon
 

On Thursday 29 January 2009 08:20, geert aarts wrote:
 Simon, thanks for your reply and your suggestions.

 I fitted the following glmm's

 gamm3-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=l
ist(code_tripnr=~1),family=poisson))

 Which worked OK (see summary below)

 I also fitted a model using quasipoisson, but that didn't help. I actually
 also thought that glmmPQL and gamm estimate the dispersion parameter and
 hence assumes a quasipoisson distribution, even if you specify poisson. Is
 that correct?

 Finally I tried fitting a model to less data, and sometimes gamm managed to
 converge (see summary below). So would it be possible to use the parameter
 estimates from the model fitted to less data as starting values for the
 gamm fitted to the full data set? Or do you have any other suggestions?

 Thanks.
 Cheers Geert






 gamm3-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=l
ist(code_tripnr=~1),f

 amily=poisson))



 iteration
 1

 iteration
 2

 iteration
 3

detach(Disc_age)

 summary(gamm3)

 Linear
 mixed-effects model fit by maximum likelihood

  Data: NULL

   AIC BIC logLik

NA  NA
 NA



 Random
 effects:

  Formula: ~1 | code_tripnr

 (Intercept) Residual

 StdDev:
 0.001391914 231.9744



 Variance
 function:

  Structure: fixed weights

  Formula: ~invwt

 Fixed
 effects: count ~ offset(offsetter) + poly(lon, 3) * poly(lat, 3)

 Value
 Std.Error   DF t-value p-value

 (Intercept)-1.582 11.96 2024 -0.13232174  0.8947

 poly(lon,
 3)1  -4.048   1397.33 2024 -0.00289673  0.9977

 poly(lon,
 3)2 -22.013699.71 2024 -0.03145996  0.9749

 poly(lon,
 3)3  -8.538593.87 2024 -0.01437683  0.9885

 poly(lat,
 3)1-109.624666.05 2024 -0.16458856  0.8693

 poly(lat,
 3)2-104.179381.37 2024 -0.27316977  0.7848

 poly(lat,
 3)3 -10.661221.93 2024 -0.04803585  0.9617

 poly(lon,
 3)1:poly(lat, 3)1  4290.737  61369.98 2024
 0.06991589  0.9443

 poly(lon,
 3)2:poly(lat, 3)1  1853.559  36835.63 2024
 0.05031972  0.9599

 poly(lon,
 3)3:poly(lat, 3)1  -240.521  25771.80 2024 -0.00933272  0.9926

 poly(lon,
 3)1:poly(lat, 3)2  2540.147  41378.38 2024
 0.06138826  0.9511

 poly(lon,
 3)1:poly(lat, 3)2  2540.147  41378.38 2024
 0.06138826  0.9511

 poly(lon,
 3)2:poly(lat, 3)2 -1803.911  21522.17
 2024 -0.08381643  0.9332

 poly(lon,
 3)3:poly(lat, 3)2  1040.858  16352.56 2024
 0.06365109  0.9493

 poly(lon,
 3)1:poly(lat, 3)3   632.587  12180.28 2024
 0.05193535  0.9586

 poly(lon,
 3)2:poly(lat, 3)3  -394.339  13088.72 2024 -0.03012818  0.9760

 poly(lon,
 3)3:poly(lat, 3)3  -543.502   6221.71 2024 -0.08735569  0.9304

  Correlation:

 (Intr) ply(ln,3)1
 ply(ln,3)2 ply(ln,3)3 ply(lt,3)1

 poly(lon,
 3)10.889

 poly(lon,
 3)20.938  0.878

 poly(lon,
 3)30.843  0.981
 0.792

 poly(lat,
 3)1   -0.829 -0.949 -0.906
 -0.882

 poly(lat,
 3)20.859  0.578  0.742
 0.538 -0.474

 poly(lat,
 3)3   -0.552 -0.783 -0.579
 -0.756  0.837

 poly(lon,
 3)1:poly(lat, 3)1 -0.947 -0.974
 -0.940 -0.940  0.925

 poly(lon,
 3)2:poly(lat, 3)1 -0.934 -0.950
 -0.857 -0.929  0.881

 poly(lon,
 3)3:poly(lat, 3)1 -0.818 -0.963
 -0.866 -0.945  0.931

 poly(lon,
 3)1:poly(lat, 3)2  0.808  0.975
 0.784  0.968 -0.928

 poly(lon,
 3)2:poly(lat, 3)2  0.737  0.575
 0.853  0.465 -0.659

 poly(lon,
 3)3:poly(lat, 3)2  0.735  0.896
 0.647  0.938 -0.765

 poly(lon,
 3)1:poly(lat, 3)3 -0.794 -0.592
 -0.823 -0.518  0.591

 poly(lon,
 3)2:poly(lat, 3)3 -0.542 -0.737
 -0.419 -0.781  0.635

 poly(lon,
 3)3:poly(lat, 3)3 -0.398 -0.383
 -0.534 -0.334  0.425

 ply(lt,3)2
 ply(lt,3)3 p(,3)1:(,3)1 p(,3)2:(,3)1

 poly(lon,
 3)1

 poly(lon,
 3)2

 poly(lon,
 3)3

 poly(lat,
 3)1

 poly(lat,
 3)2

 poly(lat,
 3)3   -0.136

 poly(lon,
 3)1:poly(lat, 3)1 -0.708

Re: [R] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread davidr

I apologize for posting a wrong opinion; I should of course have checked
before posting.

Henrik's examples illustrate something I had never realized before, and
it really surprised me!
Where can I read the technical details of this scoping aspect of 'for'
as it still presents
some subtle puzzles for me. For example:
 for (ii in 1:3) { print(ii); if(ii==1){ii - 20}; print(ii); print(ii
- ii + 1); print(ii)}
[1] 1
[1] 20
[1] 21
[1] 21
[1] 2
[1] 2
[1] 3
[1] 3
[1] 3
[1] 3
[1] 4
[1] 4
Why does R treat ii differently after the 'if' in the first and
subsequent iterations?
And if the loop variable does not exist before the 'for', why is it
created in the parent(?) environment at all?
(I.e, if ii did not exist before the for loop, it does and has value 5
after. Wouldn't strict
scoping mean that ii exists only within the for loop?)

I generally don't try to change the loop variable's value inside a loop,
but coming from C
or other languages, it seems natural to do so in certain circumstances.
And the examples are
certainly bad programming style. But I remain confused. I assume that
'while' loops have different scoping
rules?

Thanks,

-- David


-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On
Behalf Of Henrik Bengtsson
Sent: Wednesday, January 28, 2009 5:08 PM
To: David Reiner dav...@rhotrading.com
Cc: SnowManPaddington; r-help@r-project.org
Subject: Re: [R] - Re: for/if loop

On Wed, Jan 28, 2009 at 2:41 PM,  dav...@rhotrading.com wrote:
 Well, maybe you are just bad at typing then ;-)

 The lines rr==ii, pp==pp+1, etc. are not setting rr and pp but
comparing
 them.
 Probably you want rr - ii and pp - pp+1, etc.
 And the last line of your loop 'ii=ii+1' means that,
 since the for statement is already incrementing ii,
 you are incrementing it twice and skipping the even indices. Omit this
 line probably.

That is actually not the case (because of the scoping rules for for(),
I think).  Example:

 for (ii in 1:5) { print(ii); ii - ii + 1; }
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

Another counter intuitive (though it isn't) example:

for (ii in 1:3) {
  cat(Outer ii:,ii,\n);
  for (ii in ii:3) {
cat(  Inner ii:,ii,\n);
  }
}

Outer ii: 1
  Inner ii: 1
  Inner ii: 2
  Inner ii: 3
Outer ii: 2
  Inner ii: 2
  Inner ii: 3
Outer ii: 3
  Inner ii: 3

My $.02

/Henrik

 You are also forgetting(?) the operator precedence in
 sum(lselb1[rr:ii-1]) and similar lines.
 Note that this is equivalent to sum(lselb1[(rr-1):(ii-1)]); is that
what
 you meant?
 Or did you want sum(lselb1[rr:(ii-1)])?
 You are changing some variables but not asking R to print anything as
 far as I can see.
 To see the results, ask R to print hll.

 HTH,
 -- David

 -Original Message-
 From: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org]
 On Behalf Of SnowManPaddington
 Sent: Wednesday, January 28, 2009 3:59 PM
 To: r-help@r-project.org
 Subject: - Re: [R] for/if loop


 Hi ya, I've revised the code (and finally know what I m doing.. :-D)

 The good news is.. I dont get any error message, but the bad news is
the
 following optim generate no results. I still think there is something
to
 do
 with my loop... can anyone advice? Thanks again!!!



 pp=1
 rr=1

 for (ii in 1:n){
if (!(panel[ii] == pp)){
hll[pp,1] == sum(lselb1[rr:ii-1])
hll[pp,2] == sum(lselb2[rr:ii-1])
rr==ii
pp==pp+1
}

if (ii==n){
hll[pp,1] == sum(lselb1[rr:ii])
hll[pp,2] == sum(lselb2[rr:ii])
rr==ii
pp==pp+1
}
ii=ii+1
 }





 pp=1
 rr=1

 for (ii in 1:n){
if (!(panel[ii] == pp)){
hll[pp,1] == sum(lselb1[rr:ii-1])
hll[pp,2] == sum(lselb2[rr:ii-1])
rr==ii
pp==pp+1
}

if (ii==n){
hll[pp,1] == sum(lselb1[rr:ii])
hll[pp,2] == sum(lselb2[rr:ii])
rr==ii
pp==pp+1
}
ii=ii+1
 }





 SnowManPaddington wrote:

 Hi, it's my first time to write a loop with R for my homework. This
 loop
 is part of the function. I wanna assign values for hll according to
 panel
 [ii,1]=pp. I didn't get any error message in this part. but then when
 I
 further calculate another stuff with hll, the function can't return.
I
 think it must be some problem in my loop. Probably something stupid
or
 easy. But I tried to look for previous posts in forum and read R
 language
 help. But none can help.. Thanks!



 for (ii in 1:100){
   for (pp in 1:pp+1){
   for (rr in 1:rr+1){
   if (panel[ii,1]!=pp)
   {
   hll(pp,1)=ColSums(lselb1(rr:ii-1,1))
   hll(pp,2)=ColSums(lselb2(rr:ii-1,1))
   rr=ii
   pp=pp+1
   }
   else

Re: [R] Inconsistency in F values from dropterm and anova

2009-01-29 Thread Prof Brian Ripley

Perhaps if you followed the posting guide and did not send HTML mail 
your tables would be readable.  But your anova call seems wrong, as 
you have the models in decreasing not increasing order.


The correct result is (computing F test by hand)


23.1814/93

[1] 0.2492624

1.5475/3

[1] 0.5158333

0.5158333/0.2492624

[1] 2.069439

so the incorrect use of anova has given you an invalid result 
(although I don't see immediately what it has done)


On Thu, 29 Jan 2009, IMS wrote:



Hi,

I'm working on fitting a glm model to my data using Gamma error structure
and reciprocal link.  I've been using dropterm (MASS) in the model
simplification process, but the F values from analysis of deviance tables
reported by dropterm and anova functions are different - sometimes
significantly so.  However, the reported residual deviances, degrees of
freedom, etc. are not different.  I don't know how to calculate F values
from deviance tables (and can't seem to find anything suggesting how), and
so haven't been able to calculate F for myself to see which is more
accurate.  Below is an example of the code, and I'm using R version 2.8.0.


model1=glm(diff2^(0.491)~mtype*morder,family=Gamma)
dropterm(model1,test=F)

Single term deletions

Model:
diff2^(0.491) ~ mtype * morder
   Df  Deviance AIC   F value   Pr(F)
  23.181   -16.813
mtype:morder  3   24.729   -13.741  2.0694   0.1096


model2=update(model1,~.-mtype:morder)
anova(model1,model2,test=F)

Analysis of Deviance Table

Model 1: diff2^(0.491) ~ mtype * morder
Model 2: diff2^(0.491) ~ mtype + morder
 Resid. Df Resid. Dev Df Deviance  F  Pr(F)
19323.1814
29624.7288 -3  -1.5475 3.0241 0.03352 *
---
Signif. codes:  0 ???***??? 0.001 ???**??? 0.01 ???*??? 0.05 ???.??? 0.1 ??? 
??? 1

--
View this message in context: 
http://www.nabble.com/Inconsistency-in-F-values-from-dropterm-and-anova-tp21725668p21725668.html
Sent from the R help mailing list archive at Nabble.com.

[[alternative HTML version deleted]]




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

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

Re: [R] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread Duncan Murdoch


On 1/29/2009 10:39 AM, dav...@rhotrading.com wrote:

I apologize for posting a wrong opinion; I should of course have checked
before posting.

Henrik's examples illustrate something I had never realized before, and
it really surprised me!
Where can I read the technical details of this scoping aspect of 'for'
as it still presents
some subtle puzzles for me. For example:

for (ii in 1:3) { print(ii); if(ii==1){ii - 20}; print(ii); print(ii

- ii + 1); print(ii)}
[1] 1
[1] 20
[1] 21
[1] 21
[1] 2
[1] 2
[1] 3
[1] 3
[1] 3
[1] 3
[1] 4
[1] 4
Why does R treat ii differently after the 'if' in the first and
subsequent iterations?


I don't see that it does.  First time through, it prints 1, 20, 21, 21. 
 Second time through (without ii==1 being true), it prints 2, 2, 3, 3.

Third time through it prints 3,3,4,4.  What's different?



And if the loop variable does not exist before the 'for', why is it
created in the parent(?) environment at all?


It's created in the current evaluation frame, because that's where 
everything gets created unless you work hard to have it created 
somewhere else.



(I.e, if ii did not exist before the for loop, it does and has value 5
after. Wouldn't strict
scoping mean that ii exists only within the for loop?)


R doesn't have scoping as strict as that.  For loops don't create their 
own evaluation frames.  If they did, you could not do something like this:


 x - 0
 for (i in 1:10) {
  x - x + i
 }

as a slow way to do sum(1:10), because the assignment within the loop 
would take place in a local evaluation frame, and be lost afterwards.



I generally don't try to change the loop variable's value inside a loop,
but coming from C
or other languages, it seems natural to do so in certain circumstances.


R is not C.  Feel free to do what you like to the variable within the 
loop (though you might cause Luke to grind his teeth when it messes up 
his compiler).  It will be set to the next value in the 1:10 vector the 
next time it comes back to the top.



And the examples are
certainly bad programming style. But I remain confused. I assume that
'while' loops have different scoping
rules?


No, the scoping rules in R are quite simple and consistent.  for and 
while have the same rules.


Duncan Murdoch



Thanks,

-- David


-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On
Behalf Of Henrik Bengtsson
Sent: Wednesday, January 28, 2009 5:08 PM
To: David Reiner dav...@rhotrading.com
Cc: SnowManPaddington; r-help@r-project.org
Subject: Re: [R] - Re: for/if loop

On Wed, Jan 28, 2009 at 2:41 PM,  dav...@rhotrading.com wrote:

Well, maybe you are just bad at typing then ;-)

The lines rr==ii, pp==pp+1, etc. are not setting rr and pp but

comparing

them.
Probably you want rr - ii and pp - pp+1, etc.
And the last line of your loop 'ii=ii+1' means that,
since the for statement is already incrementing ii,
you are incrementing it twice and skipping the even indices. Omit this
line probably.


That is actually not the case (because of the scoping rules for for(),
I think).  Example:


for (ii in 1:5) { print(ii); ii - ii + 1; }

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

Another counter intuitive (though it isn't) example:

for (ii in 1:3) {
  cat(Outer ii:,ii,\n);
  for (ii in ii:3) {
cat(  Inner ii:,ii,\n);
  }
}

Outer ii: 1
  Inner ii: 1
  Inner ii: 2
  Inner ii: 3
Outer ii: 2
  Inner ii: 2
  Inner ii: 3
Outer ii: 3
  Inner ii: 3

My $.02

/Henrik


You are also forgetting(?) the operator precedence in
sum(lselb1[rr:ii-1]) and similar lines.
Note that this is equivalent to sum(lselb1[(rr-1):(ii-1)]); is that

what

you meant?
Or did you want sum(lselb1[rr:(ii-1)])?
You are changing some variables but not asking R to print anything as
far as I can see.
To see the results, ask R to print hll.

HTH,
-- David

-Original Message-
From: r-help-boun...@r-project.org

[mailto:r-help-boun...@r-project.org]

On Behalf Of SnowManPaddington
Sent: Wednesday, January 28, 2009 3:59 PM
To: r-help@r-project.org
Subject: - Re: [R] for/if loop


Hi ya, I've revised the code (and finally know what I m doing.. :-D)

The good news is.. I dont get any error message, but the bad news is

the

following optim generate no results. I still think there is something

to

do
with my loop... can anyone advice? Thanks again!!!



pp=1
rr=1

for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){
   hll[pp,1] == sum(lselb1[rr:ii])
   hll[pp,2] == sum(lselb2[rr:ii])
   rr==ii
   pp==pp+1
   }
   ii=ii+1
}





pp=1
rr=1

for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){

Re: [R] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread Patrick Burns


Certainly not a complete description, but
'The R Inferno' talks about this on page 62.


Patrick Burns
patr...@burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of The R Inferno and A Guide for the Unwilling S User)

dav...@rhotrading.com wrote:

I apologize for posting a wrong opinion; I should of course have checked
before posting.

Henrik's examples illustrate something I had never realized before, and
it really surprised me!
Where can I read the technical details of this scoping aspect of 'for'
as it still presents
some subtle puzzles for me. For example:
  

for (ii in 1:3) { print(ii); if(ii==1){ii - 20}; print(ii); print(ii


- ii + 1); print(ii)}
[1] 1
[1] 20
[1] 21
[1] 21
[1] 2
[1] 2
[1] 3
[1] 3
[1] 3
[1] 3
[1] 4
[1] 4
Why does R treat ii differently after the 'if' in the first and
subsequent iterations?
And if the loop variable does not exist before the 'for', why is it
created in the parent(?) environment at all?
(I.e, if ii did not exist before the for loop, it does and has value 5
after. Wouldn't strict
scoping mean that ii exists only within the for loop?)

I generally don't try to change the loop variable's value inside a loop,
but coming from C
or other languages, it seems natural to do so in certain circumstances.
And the examples are
certainly bad programming style. But I remain confused. I assume that
'while' loops have different scoping
rules?

Thanks,

-- David


-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On
Behalf Of Henrik Bengtsson
Sent: Wednesday, January 28, 2009 5:08 PM
To: David Reiner dav...@rhotrading.com
Cc: SnowManPaddington; r-help@r-project.org
Subject: Re: [R] - Re: for/if loop

On Wed, Jan 28, 2009 at 2:41 PM,  dav...@rhotrading.com wrote:
  

Well, maybe you are just bad at typing then ;-)

The lines rr==ii, pp==pp+1, etc. are not setting rr and pp but


comparing
  

them.
Probably you want rr - ii and pp - pp+1, etc.
And the last line of your loop 'ii=ii+1' means that,
since the for statement is already incrementing ii,
you are incrementing it twice and skipping the even indices. Omit this
line probably.



That is actually not the case (because of the scoping rules for for(),
I think).  Example:

  

for (ii in 1:5) { print(ii); ii - ii + 1; }


[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

Another counter intuitive (though it isn't) example:

for (ii in 1:3) {
  cat(Outer ii:,ii,\n);
  for (ii in ii:3) {
cat(  Inner ii:,ii,\n);
  }
}

Outer ii: 1
  Inner ii: 1
  Inner ii: 2
  Inner ii: 3
Outer ii: 2
  Inner ii: 2
  Inner ii: 3
Outer ii: 3
  Inner ii: 3

My $.02

/Henrik

  

You are also forgetting(?) the operator precedence in
sum(lselb1[rr:ii-1]) and similar lines.
Note that this is equivalent to sum(lselb1[(rr-1):(ii-1)]); is that


what
  

you meant?
Or did you want sum(lselb1[rr:(ii-1)])?
You are changing some variables but not asking R to print anything as
far as I can see.
To see the results, ask R to print hll.

HTH,
-- David

-Original Message-
From: r-help-boun...@r-project.org


[mailto:r-help-boun...@r-project.org]
  

On Behalf Of SnowManPaddington
Sent: Wednesday, January 28, 2009 3:59 PM
To: r-help@r-project.org
Subject: - Re: [R] for/if loop


Hi ya, I've revised the code (and finally know what I m doing.. :-D)

The good news is.. I dont get any error message, but the bad news is


the
  

following optim generate no results. I still think there is something


to
  

do
with my loop... can anyone advice? Thanks again!!!



pp=1
rr=1

for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){
   hll[pp,1] == sum(lselb1[rr:ii])
   hll[pp,2] == sum(lselb2[rr:ii])
   rr==ii
   pp==pp+1
   }
   ii=ii+1
}





pp=1
rr=1

for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){
   hll[pp,1] == sum(lselb1[rr:ii])
   hll[pp,2] == sum(lselb2[rr:ii])
   rr==ii
   pp==pp+1
   }
   ii=ii+1
}





SnowManPaddington wrote:


Hi, it's my first time to write a loop with R for my homework. This
  

loop


is part of the function. I wanna assign values for hll according to
  

panel


[ii,1]=pp. I didn't get any error message in this part. but then when
  

I


further calculate another stuff with hll, the function can't return.
  

I
  

think it must be some problem in my loop. Probably something stupid
  

or
  

easy. But I tried to look for previous posts in forum and read R
  

language


help. But none can help.. Thanks!



for (ii in 1:100){
  for (pp

[R] Help

2009-01-29 Thread Alexandra Almeida

Hi everybody!

I´m with a problem that probably is easy for you, but I really don´t know
how to solve.
On the following script:

for(j in 1:length(limiares))
{
   excessos-limiares[j]-estacao[estacaolimiares[j]]
   par.ests-gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.est

GOF.test-ks.test(excessos,pgpd,xi=par.ests[1],beta=par.ests[2])$p.value

tabs[j,]-c(par.ests,gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.ses,GOF.test,length(excessos))
}

I´ve found the error for some values of i: Erro em ks.test(excessos, pgpd,
xi = par.ests[1], beta = par.ests[2]) : NA/NaN/Inf em chamada de função
externa (argumento 1)

*My question is: This warning stop the for and I don´t want it, is there
some way to continue the for, and for the cases where the function cannot
calculate the ks.test for the i just leave a NA as answer???*

Thank you *very much*!!!

Alexandra Almeida



-- 
Alexandra R M de Almeida

[[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] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread Duncan Murdoch


On 1/29/2009 11:06 AM, Patrick Burns wrote:

Certainly not a complete description, but
'The R Inferno' talks about this on page 62.


Hmmm, I don't think I agree with that description.  There's only one 
object named i in that example (which was


for (i in 1:6) {
  cat('\n i is', i, '\n\n')
  for (i in 1:4) cat('i is', i,'\n')
}

for the lazy).  It gets set to 1, then to each of 1:4, then to 2, then 
to 1:4, etc.  If you print the value after the inner loop, it will 
always be 4, you won't go back to the supposed other i.  For example,


 for (i in 1:6) {
+   cat('\n i is', i, '\n\n')
+   for (i in 1:4) cat('i is', i,'\n')
+   cat('outer i is', i, '\n')
+ }

 i is 1

i is 1
i is 2
i is 3
i is 4
outer i is 4

 i is 2

i is 1
i is 2
i is 3
i is 4
outer i is 4

 i is 3

i is 1
i is 2
i is 3
i is 4
outer i is 4

Duncan Murdoch





Patrick Burns
patr...@burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of The R Inferno and A Guide for the Unwilling S User)

dav...@rhotrading.com wrote:

I apologize for posting a wrong opinion; I should of course have checked
before posting.

Henrik's examples illustrate something I had never realized before, and
it really surprised me!
Where can I read the technical details of this scoping aspect of 'for'
as it still presents
some subtle puzzles for me. For example:
  

for (ii in 1:3) { print(ii); if(ii==1){ii - 20}; print(ii); print(ii


- ii + 1); print(ii)}
[1] 1
[1] 20
[1] 21
[1] 21
[1] 2
[1] 2
[1] 3
[1] 3
[1] 3
[1] 3
[1] 4
[1] 4
Why does R treat ii differently after the 'if' in the first and
subsequent iterations?
And if the loop variable does not exist before the 'for', why is it
created in the parent(?) environment at all?
(I.e, if ii did not exist before the for loop, it does and has value 5
after. Wouldn't strict
scoping mean that ii exists only within the for loop?)

I generally don't try to change the loop variable's value inside a loop,
but coming from C
or other languages, it seems natural to do so in certain circumstances.
And the examples are
certainly bad programming style. But I remain confused. I assume that
'while' loops have different scoping
rules?

Thanks,

-- David


-Original Message-
From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On
Behalf Of Henrik Bengtsson
Sent: Wednesday, January 28, 2009 5:08 PM
To: David Reiner dav...@rhotrading.com
Cc: SnowManPaddington; r-help@r-project.org
Subject: Re: [R] - Re: for/if loop

On Wed, Jan 28, 2009 at 2:41 PM,  dav...@rhotrading.com wrote:
  

Well, maybe you are just bad at typing then ;-)

The lines rr==ii, pp==pp+1, etc. are not setting rr and pp but


comparing
  

them.
Probably you want rr - ii and pp - pp+1, etc.
And the last line of your loop 'ii=ii+1' means that,
since the for statement is already incrementing ii,
you are incrementing it twice and skipping the even indices. Omit this
line probably.



That is actually not the case (because of the scoping rules for for(),
I think).  Example:

  

for (ii in 1:5) { print(ii); ii - ii + 1; }


[1] 1
[1] 2
[1] 3
[1] 4
[1] 5

Another counter intuitive (though it isn't) example:

for (ii in 1:3) {
  cat(Outer ii:,ii,\n);
  for (ii in ii:3) {
cat(  Inner ii:,ii,\n);
  }
}

Outer ii: 1
  Inner ii: 1
  Inner ii: 2
  Inner ii: 3
Outer ii: 2
  Inner ii: 2
  Inner ii: 3
Outer ii: 3
  Inner ii: 3

My $.02

/Henrik

  

You are also forgetting(?) the operator precedence in
sum(lselb1[rr:ii-1]) and similar lines.
Note that this is equivalent to sum(lselb1[(rr-1):(ii-1)]); is that


what
  

you meant?
Or did you want sum(lselb1[rr:(ii-1)])?
You are changing some variables but not asking R to print anything as
far as I can see.
To see the results, ask R to print hll.

HTH,
-- David

-Original Message-
From: r-help-boun...@r-project.org


[mailto:r-help-boun...@r-project.org]
  

On Behalf Of SnowManPaddington
Sent: Wednesday, January 28, 2009 3:59 PM
To: r-help@r-project.org
Subject: - Re: [R] for/if loop


Hi ya, I've revised the code (and finally know what I m doing.. :-D)

The good news is.. I dont get any error message, but the bad news is


the
  

following optim generate no results. I still think there is something


to
  

do
with my loop... can anyone advice? Thanks again!!!



pp=1
rr=1

for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){
   hll[pp,1] == sum(lselb1[rr:ii])
   hll[pp,2] == sum(lselb2[rr:ii])
   rr==ii
   pp==pp+1
   }
   ii=ii+1
}





pp=1
rr=1

for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){
   hll[pp,1] ==

Re: [R] Help

2009-01-29 Thread patricia garcía gonzález


Hi,


Test - function(x, y, paired) {
   res - try( ks.test(...), TRUE )
  if( inherits( res, try-error ) )
  res - NA
  res
}


Regards


patricia


 Date: Thu, 29 Jan 2009 13:22:11 -0300
 From: alexandra...@gmail.com
 To: r-help@r-project.org
 Subject: [R] Help
 
 Hi everybody!
 
 I´m with a problem that probably is easy for you, but I really don´t know
 how to solve.
 On the following script:
 
 for(j in 1:length(limiares))
 {
excessos-limiares[j]-estacao[estacaolimiares[j]]
par.ests-gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.est
 
 GOF.test-ks.test(excessos,pgpd,xi=par.ests[1],beta=par.ests[2])$p.value
 
 tabs[j,]-c(par.ests,gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.ses,GOF.test,length(excessos))
 }
 
 I´ve found the error for some values of i: Erro em ks.test(excessos, pgpd,
 xi = par.ests[1], beta = par.ests[2]) : NA/NaN/Inf em chamada de função
 externa (argumento 1)
 
 *My question is: This warning stop the for and I don´t want it, is there
 some way to continue the for, and for the cases where the function cannot
 calculate the ks.test for the i just leave a NA as answer???*
 
 Thank you *very much*!!!
 
 Alexandra Almeida
 
 
 
 -- 
 Alexandra R M de Almeida
 
   [[alternative HTML version deleted]]
 

_


[[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] Standard Errors of Coefficients

2009-01-29 Thread Andreas Klein

Hello.

I estimated a VAR(1) TSmodel (var_1) with estMaxLik from the dse1 package given 
a TSdata object (mydata).

est.model - estMaxLik(var_1,mydata)

How can I obtain the standard errors for the four coefficients of the estimated 
model to check for significance? - Is it yet calculated and I can obtain it by 
est.model$  ?


Thank you in advance.

Regards,
Andreas.




__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] bootstrapping in regression

2009-01-29 Thread Thomas Mang


Hi,

Please apologize if my questions sounds somewhat 'stupid' to the trained 
and experienced statisticians of you. Also I am not sure if I used all 
terms correctly, if not then corrections are welcome.


I have asked myself the following question regarding bootstrapping in 
regression:
Say for whatever reason one does not want to take the p-values for 
regression coefficients from the established test statistics 
distributions (t-distr for individual coefficients, F-values for 
whole-model-comparisons), but instead apply a more robust approach by 
bootstrapping.


In the simple linear regression case, one possibility is to randomly 
rearrange the X/Y data pairs, estimate the model and take the 
beta1-coefficient. Do this many many times, and so derive the null 
distribution for beta1. Finally compare beta1 for the observed data 
against this null-distribution.


What I now wonder is how the situation looks like in the multiple 
regression case. Assume there are two predictors, X1 and X2. Is it then 
possible to do the same, but just only rearranging the values of one 
predictor (the one of interest) at a time? Say I want again to test 
beta1. Is it then valid to many times randomly rearrange the X1 data 
(and keeping Y and X2 as observed), fit the model, take the beta1 
coefficient, and finally compare the beta1 of the observed data against 
the distributions of these beta1s ?
For X2, do the same, randomly rearrange X2 all the time while keeping Y 
and X1 as observed etc.

Is this valid ?

Second, if this is valid for the 'normal', fixed-effects only 
regression, is it also valid to derive null distributions for the 
regression coefficients of the fixed effects in a mixed model this way? 
Or does the quite different parameters estimation calculation forbid 
this approach (Forbid in the sense of bogus outcome) ?


Thanks, Thomas

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

Re: [R] Help

2009-01-29 Thread rkevinburton

You could look into ' try' and set it up to catch errors and do the appropriate 
thing in your error handler. I don't have the exact syntax at hand right now 
but looking at ?try or ?tryCatch I think will do what you want.

Kevin

 Alexandra Almeida alexandra...@gmail.com wrote: 
 Hi everybody!
 
 I´m with a problem that probably is easy for you, but I really don´t know
 how to solve.
 On the following script:
 
 for(j in 1:length(limiares))
 {
excessos-limiares[j]-estacao[estacaolimiares[j]]
par.ests-gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.est
 
 GOF.test-ks.test(excessos,pgpd,xi=par.ests[1],beta=par.ests[2])$p.value
 
 tabs[j,]-c(par.ests,gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.ses,GOF.test,length(excessos))
 }
 
 I´ve found the error for some values of i: Erro em ks.test(excessos, pgpd,
 xi = par.ests[1], beta = par.ests[2]) : NA/NaN/Inf em chamada de função
 externa (argumento 1)
 
 *My question is: This warning stop the for and I don´t want it, is there
 some way to continue the for, and for the cases where the function cannot
 calculate the ks.test for the i just leave a NA as answer???*
 
 Thank you *very much*!!!
 
 Alexandra Almeida
 
 
 
 -- 
 Alexandra R M de Almeida
 
   [[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] Taking the min of each row in a matrix

2009-01-29 Thread Will Glass-Husain

Hi,

I'm a new user. I've been reading through the manual and looking at
various examples but am still trying to make sense of the most
efficient ways to handle matrices of data.

If I have a 2D matrix of data, how do I get the mean, min, max value
of each row?  I see the mean function on a matrix will give me
averages by row, but min and max give me the value for the entire
matrix.  I want the min (for example) of each row.  pmin looks useful,
but I can't seem to get the syntax right to apply to each column.

Right now I am doing this.  Is there a one-liner that would work instead?

minResult - vector (mode=list, length=rowCount)
for (row in 1:rowCount)
{
minResult[[row]] - min(corResult[row], na.rm = TRUE)   
}

Thanks in advance.

WILL

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread Henrik Bengtsson

Thanks for straighten this out.  Sorry for my misleading suggestion
that special scoping rules comes into play; it is just about
reassignments at the beginning of each loop.  Here is a slightly
better illustration:

ii - start;
cat(ii:,ii,\n);
for (ii in 1:2) {
 cat(Outer ii:,ii,\n);
 for (ii in c(a, b, c)) {
   cat(  Inner ii:,ii,\n);
 }
 cat(At outer ii:,ii,\n);
}
cat(ii:,ii,\n);

ii: start
Outer ii: 1
  Inner ii: a
  Inner ii: b
  Inner ii: c
At outer ii: c
Outer ii: 2
  Inner ii: a
  Inner ii: b
  Inner ii: c
At outer ii: c
ii: c

/Henrik

PS. About the double-letter index (e.g. ii vs. i); A few years ago
someone suggested me to use this, because it is much easier to search
for 'ii' in an editor compared with a single-letter 'i'.  So true.  I
made the move immediately.

On Thu, Jan 29, 2009 at 8:06 AM, Patrick Burns pbu...@pburns.seanet.com wrote:
 Certainly not a complete description, but
 'The R Inferno' talks about this on page 62.


 Patrick Burns
 patr...@burns-stat.com
 +44 (0)20 8525 0696
 http://www.burns-stat.com
 (home of The R Inferno and A Guide for the Unwilling S User)

 dav...@rhotrading.com wrote:

 I apologize for posting a wrong opinion; I should of course have checked
 before posting.

 Henrik's examples illustrate something I had never realized before, and
 it really surprised me!
 Where can I read the technical details of this scoping aspect of 'for'
 as it still presents
 some subtle puzzles for me. For example:


 for (ii in 1:3) { print(ii); if(ii==1){ii - 20}; print(ii); print(ii


 - ii + 1); print(ii)}
 [1] 1
 [1] 20
 [1] 21
 [1] 21
 [1] 2
 [1] 2
 [1] 3
 [1] 3
 [1] 3
 [1] 3
 [1] 4
 [1] 4
 Why does R treat ii differently after the 'if' in the first and
 subsequent iterations?
 And if the loop variable does not exist before the 'for', why is it
 created in the parent(?) environment at all?
 (I.e, if ii did not exist before the for loop, it does and has value 5
 after. Wouldn't strict
 scoping mean that ii exists only within the for loop?)

 I generally don't try to change the loop variable's value inside a loop,
 but coming from C
 or other languages, it seems natural to do so in certain circumstances.
 And the examples are
 certainly bad programming style. But I remain confused. I assume that
 'while' loops have different scoping
 rules?

 Thanks,

 -- David


 -Original Message-
 From: henrik.bengts...@gmail.com [mailto:henrik.bengts...@gmail.com] On
 Behalf Of Henrik Bengtsson
 Sent: Wednesday, January 28, 2009 5:08 PM
 To: David Reiner dav...@rhotrading.com
 Cc: SnowManPaddington; r-help@r-project.org
 Subject: Re: [R] - Re: for/if loop

 On Wed, Jan 28, 2009 at 2:41 PM,  dav...@rhotrading.com wrote:


 Well, maybe you are just bad at typing then ;-)

 The lines rr==ii, pp==pp+1, etc. are not setting rr and pp but


 comparing


 them.
 Probably you want rr - ii and pp - pp+1, etc.
 And the last line of your loop 'ii=ii+1' means that,
 since the for statement is already incrementing ii,
 you are incrementing it twice and skipping the even indices. Omit this
 line probably.


 That is actually not the case (because of the scoping rules for for(),
 I think).  Example:



 for (ii in 1:5) { print(ii); ii - ii + 1; }


 [1] 1
 [1] 2
 [1] 3
 [1] 4
 [1] 5

 Another counter intuitive (though it isn't) example:

 for (ii in 1:3) {
  cat(Outer ii:,ii,\n);
  for (ii in ii:3) {
cat(  Inner ii:,ii,\n);
  }
 }

 Outer ii: 1
  Inner ii: 1
  Inner ii: 2
  Inner ii: 3
 Outer ii: 2
  Inner ii: 2
  Inner ii: 3
 Outer ii: 3
  Inner ii: 3

 My $.02

 /Henrik



 You are also forgetting(?) the operator precedence in
 sum(lselb1[rr:ii-1]) and similar lines.
 Note that this is equivalent to sum(lselb1[(rr-1):(ii-1)]); is that


 what


 you meant?
 Or did you want sum(lselb1[rr:(ii-1)])?
 You are changing some variables but not asking R to print anything as
 far as I can see.
 To see the results, ask R to print hll.

 HTH,
 -- David

 -Original Message-
 From: r-help-boun...@r-project.org


 [mailto:r-help-boun...@r-project.org]


 On Behalf Of SnowManPaddington
 Sent: Wednesday, January 28, 2009 3:59 PM
 To: r-help@r-project.org
 Subject: - Re: [R] for/if loop


 Hi ya, I've revised the code (and finally know what I m doing.. :-D)

 The good news is.. I dont get any error message, but the bad news is


 the


 following optim generate no results. I still think there is something


 to


 do
 with my loop... can anyone advice? Thanks again!!!



 pp=1
 rr=1

 for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])
   hll[pp,2] == sum(lselb2[rr:ii-1])
   rr==ii
   pp==pp+1
   }

   if (ii==n){
   hll[pp,1] == sum(lselb1[rr:ii])
   hll[pp,2] == sum(lselb2[rr:ii])
   rr==ii
   pp==pp+1
   }
   ii=ii+1
 }





 pp=1
 rr=1

 for (ii in 1:n){
   if (!(panel[ii] == pp)){
   hll[pp,1] == sum(lselb1[rr:ii-1])

Re: [R] Dynamic random effects model

2009-01-29 Thread Joseph Magagnoli

I am trying to estimate a model of third party intervention into civil war.
The data is panel data structure.  The unit of analysis
is civil war year.  For each civil war year my dependent variable is coded
0=no intervention and 1=intervention.  I want to use a
lagged dependent variable as an independent variable and i am including
other variables such as type of war conventional=0,1 dummy
irregular 0,1 dummy,,, other dummy variables such as cold war period, and
other variables such as state strength .   Some of my independent variables
are time invariant or slow moving.   Because i want to include lagged
dependent variable rules out a fixed effect,
therefore I was thinking of using random effects.Any suggestion on how
to model this?

I am inexperienced with these models, I will appreciate any help I can get

Thank you

jcm

On Thu, Jan 29, 2009 at 8:04 AM, Ben Bolker bol...@ufl.edu wrote:

  Joseph Magagnoli jcm331 at gmail.com writes:

 
  All R experts,
  How do I fit a dynamic Random effects model with a binary dependent
 variable
  in R
  Thanks
  JCM
 

  You haven't given us nearly enough information to go on.
 If you're talking about something like a state-space model with
 a binary response, I would probably say your best bet is
 a Bayesian approach, prob. via JAGS/R2jags or WinBUGS/R2WinBUGS.
 (A lot) more context will give a higher probability of a useful
 answer.

  good luck
Ben Bolker

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.htmlhttp://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] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread Wacek Kusnierczyk

Duncan Murdoch wrote:
 On 1/29/2009 10:39 AM, dav...@rhotrading.com wrote:

snip


 And if the loop variable does not exist before the 'for', why is it
 created in the parent(?) environment at all?

 It's created in the current evaluation frame, because that's where
 everything gets created unless you work hard to have it created
 somewhere else.

... and it's by no means idiosyncratic to r; many dynamic languages do
it this way.


 (I.e, if ii did not exist before the for loop, it does and has value 5
 after. Wouldn't strict
 scoping mean that ii exists only within the for loop?)

 R doesn't have scoping as strict as that.  For loops don't create
 their own evaluation frames.  If they did, you could not do something
 like this:

  x - 0
  for (i in 1:10) {
   x - x + i
  }

 as a slow way to do sum(1:10), because the assignment within the loop
 would take place in a local evaluation frame, and be lost afterwards.

 I generally don't try to change the loop variable's value inside a loop,
 but coming from C
 or other languages, it seems natural to do so in certain circumstances.

 R is not C.  Feel free to do what you like to the variable within the
 loop (though you might cause Luke to grind his teeth when it messes up
 his compiler).  It will be set to the next value in the 1:10 vector
 the next time it comes back to the top.

in the iso c99 standard the following is illegal:

for (int i = 0; i  10; i++) ...

and you have to declare the iterator variable in advance:

int i;
for (i = 0; i  10; i++) ...

which effectively works as the r code above, so no surprises if you're
coming from c99 c.  it is different in c++:

int i = 0;
for (int i = 0; i  10; i++) ...
// i still 0 here

(and this won't go in c99, again).

vQ

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Dynamic random effects model

2009-01-29 Thread Ben Bolker


  I am forwarding back to the R list in hopes of getting you a better
answer from someone else.  This sounds like it might be logistic
autoregression?

RSiteSearch(logistic autoregress*) came up with
http://finzi.psych.upenn.edu/R/library/repeated/html/gar.html

  Ben Bolker


Joseph Magagnoli wrote:
 I am trying to estimate a model of third party intervention into civil war.  
 The data is panel data structure.  The unit of analysis
 is civil war year.  For each civil war year my dependent variable is coded 
 0=no intervention and 1=intervention.  I want to use a
 lagged dependent variable as an independent variable and i am including other 
 variables such as type of war conventional=0,1 dummy
 irregular 0,1 dummy,,, other dummy variables such as cold war period, and 
 other variables such as state strength .   Some of my independent variables 
 are time invariant or slow moving.   Because i want to include lagged 
 dependent variable rules out a fixed effect,
 therefore I was thinking of using random effects.Any suggestion on how to 
 model this?
 
 I am inexperienced with these models, I will appreciate any help I can get
 
 Thank you
 
 jcm
 On Thu, Jan 29, 2009 at 8:04 AM, Ben Bolker 
 bol...@ufl.edumailto:bol...@ufl.edu wrote:
 Joseph Magagnoli jcm331 at gmail.comhttp://gmail.com/ writes:
 
 All R experts,
 How do I fit a dynamic Random effects model with a binary dependent variable
 in R
 Thanks
 JCM

 
  You haven't given us nearly enough information to go on.
 If you're talking about something like a state-space model with
 a binary response, I would probably say your best bet is
 a Bayesian approach, prob. via JAGS/R2jags or WinBUGS/R2WinBUGS.
 (A lot) more context will give a higher probability of a useful
 answer.
 
  good luck
Ben Bolker
 
 __
 R-help@r-project.orgmailto:R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
 http://www.R-project.org/posting-guide.htmlhttp://www.r-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 


-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bol...@ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc



-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bol...@ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc


signature.asc
Description: PGP signature


signature.asc
Description: OpenPGP digital signature
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Taking the min of each row in a matrix

2009-01-29 Thread anna freni sterrantino

Hello Will,

z=matrix(c(1,2,3,4,5,6), 2,byrow=T)
min=apply(z,1,min)

what do you need isapply
?apply and its family are very helpful.

apply( matrix, columns, function)
apply( matrix, rows, function)
Check   ?apply


 

Cheers

Anna



Anna Freni Sterrantino
Ph.D Student 
Department of Statistics
University of Bologna, Italy
via Belle Arti 41, 40124 BO.





Da: Will Glass-Husain wglasshus...@gmail.com
A: r-h...@stat.math.ethz.ch
Inviato: GiovedÃ¬ 29 gennaio 2009, 17:58:50
Oggetto: [R] Taking the min of each row in a matrix

Hi,

I'm a new user. I've been reading through the manual and looking at
various examples but am still trying to make sense of the most
efficient ways to handle matrices of data.

If I have a 2D matrix of data, how do I get the mean, min, max value
of each row?  I see the mean function on a matrix will give me
averages by row, but min and max give me the value for the entire
matrix.  I want the min (for example) of each row.  pmin looks useful,
but I can't seem to get the syntax right to apply to each column.

Right now I am doing this.  Is there a one-liner that would work instead?

minResult - vector (mode=list, length=rowCount)
for (row in 1:rowCount)
{
minResult[[row]] - min(corResult[row], na.rm = TRUE)
}

Thanks in advance.

WILL

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] bootstrapping in regression

2009-01-29 Thread Chuck Cleland

On 1/29/2009 11:43 AM, Thomas Mang wrote:
 Hi,
 
 Please apologize if my questions sounds somewhat 'stupid' to the trained
 and experienced statisticians of you. Also I am not sure if I used all
 terms correctly, if not then corrections are welcome.
 
 I have asked myself the following question regarding bootstrapping in
 regression:
 Say for whatever reason one does not want to take the p-values for
 regression coefficients from the established test statistics
 distributions (t-distr for individual coefficients, F-values for
 whole-model-comparisons), but instead apply a more robust approach by
 bootstrapping.
 
 In the simple linear regression case, one possibility is to randomly
 rearrange the X/Y data pairs, estimate the model and take the
 beta1-coefficient. Do this many many times, and so derive the null
 distribution for beta1. Finally compare beta1 for the observed data
 against this null-distribution.
 
 What I now wonder is how the situation looks like in the multiple
 regression case. Assume there are two predictors, X1 and X2. Is it then
 possible to do the same, but just only rearranging the values of one
 predictor (the one of interest) at a time? Say I want again to test
 beta1. Is it then valid to many times randomly rearrange the X1 data
 (and keeping Y and X2 as observed), fit the model, take the beta1
 coefficient, and finally compare the beta1 of the observed data against
 the distributions of these beta1s ?
 For X2, do the same, randomly rearrange X2 all the time while keeping Y
 and X1 as observed etc.
 Is this valid ?
 
 Second, if this is valid for the 'normal', fixed-effects only
 regression, is it also valid to derive null distributions for the
 regression coefficients of the fixed effects in a mixed model this way?
 Or does the quite different parameters estimation calculation forbid
 this approach (Forbid in the sense of bogus outcome) ?
 
 Thanks, Thomas

  Have a look at the following document by John Fox:

http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-bootstrapping.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.

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

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

Re: [R] Taking the min of each row in a matrix

2009-01-29 Thread jim holtman

?apply

e.g. for the 'min'

apply(yourMatrix, 1, min, na.rm=TRUE)

On Thu, Jan 29, 2009 at 11:58 AM, Will Glass-Husain
wglasshus...@gmail.com wrote:
 Hi,

 I'm a new user. I've been reading through the manual and looking at
 various examples but am still trying to make sense of the most
 efficient ways to handle matrices of data.

 If I have a 2D matrix of data, how do I get the mean, min, max value
 of each row?  I see the mean function on a matrix will give me
 averages by row, but min and max give me the value for the entire
 matrix.  I want the min (for example) of each row.  pmin looks useful,
 but I can't seem to get the syntax right to apply to each column.

 Right now I am doing this.  Is there a one-liner that would work instead?

 minResult - vector (mode=list, length=rowCount)
 for (row in 1:rowCount)
 {
minResult[[row]] - min(corResult[row], na.rm = TRUE)
 }

 Thanks in advance.

 WILL

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Help

2009-01-29 Thread Alexandra Almeida

Hi everybody!

I´m with a problem that probably is easy for you but I really don´t know how
to solve.
On the following script:

for(j in 1:length(limiares))
{
   excessos-limiares[j]-estacao[estacaolimiares[j]]
   par.ests-gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.est

GOF.test-ks.test(excessos,pgpd,xi=par.ests[1],beta=par.ests[2])$p.value

tabs[j,]-c(par.ests,gpd(-(estacao),threshold=-limiares[j],method=c(pwm))$par.ses,GOF.test,length(excessos))
}

I´ve found the error for some values of i: Erro em ks.test(excessos, pgpd,
xi = par.ests[1], beta = par.ests[2]) : NA/NaN/Inf em chamada de função
externa (argumento 1)

*My question is: This warning stop the for and I don´t want it, is there
some way to continue the for, and for the cases where the function cannot
calculate the ks.test for the i just leave a NA as answer???*

Thank you *very much*!!!

Alexandra Almeida






-- 
Alexandra R M de Almeida

[[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] Stacking data

2009-01-29 Thread David Winsemius


What does str(zz) tell you?

--
David Winsemius
On Jan 29, 2009, at 10:31 AM, Amit Patel wrote:


Hi

I have data in the format below

 AgeV1   V2  V3   V4
   23646 45190 50333 55166 56271
   26174 35535 38227 37911 41184
   27723  25691 25712 26144 26398

and would like to sort it as follows

 Age  values ind
2364645190V1
   26174   35535  V1
   2772325691 V2
  2719330949  V2

But i have 41 columns (age column  +  40 individuals)
I have the following script but an error is thrown up
can anyone help, where am i going wrong

zz - read.csv(Filename.csv,strip.white = TRUE)
zzz - cbind(zz[gl(nrow(zz), 1, 40*nrow(zz)), 1], stack(zz[, 2:41]))

ERROR MESSAGE
Error in data.frame(..., check.names = FALSE) :
  arguments imply differing number of rows: 3789080, 0





[[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] Accessing R as a Web service (UNCLASSIFIED)

2009-01-29 Thread Neiderer, Andrew (Civ, ARL/CISD)

Classification:  UNCLASSIFIED 
Caveats: NONE
 
I am new to R and somewhat to Web server programming.  
I am a Java programmer, however, and have done quite a bit with X3D
(extensible 3D) computer graphics.

A statistician that I work with is doing multidimensional scaling (MDS) and
provides some x,y,z's which I display using X3D.
Currently he is using Permap which more than gets the job done but is
written in Visual Basic (I believe).  I think it would be a big job to use
Permap as a Web service.  So I am looking at R (for other reasons as well).

My question is -
Can I access R via SOAP, e.g., running on a server.  Assuming I worded my
question correctly and someone understands my question, please provide URLs
. so I can get knowledgeable.

Thank you.

- Andrew M. Neiderer

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

[R] Problem VGAM and Predict

2009-01-29 Thread EinSchwede

Hello,

since I installed the package VGAM I have problems useing the predict for 
othere methods.
for example I have a model from glm and polr the command predict(model) I get 
the error: unable to find an inherited method for function predict, for 
signature polr.

Has perhaps anybody a solution, because Iwould need vglm and also other methods 
like tree in a loop.

Thanks a lot
Vinzent
-- 
NUR NOCH BIS 31.01.! GMX FreeDSL - Telefonanschluss + DSL 
für nur 16,37 EURO/mtl.!* http://dsl.gmx.de/?ac=OM.AD.PD003K11308T4569a

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread davidr

The lines below made me understand clearly. Maybe they are already in
some documentation,
but if not, it might help others to avoid my misunderstanding.
Thanks to all for the clarifications.
-- David


-Original Message-
From: Duncan Murdoch [mailto:murd...@stats.uwo.ca] 
Sent: Thursday, January 29, 2009 9:54 AM
To: David Reiner dav...@rhotrading.com
Cc: Henrik Bengtsson; r-help@r-project.org; SnowManPaddington
Subject: [SPAM] - Re: [R] scoping rules for 'for' [was - Re: for/if loop
] - Found word(s) list error in the Text body

snip
  Feel free to do what you like to the variable within the 
loop (though you might cause Luke to grind his teeth when it messes up 
his compiler).  It will be set to the next value in the 1:10 vector the 
next time it comes back to the top.

snip

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.

[R] Unicode

2009-01-29 Thread R Heberto Ghezzo, Dr

Hello, can somebody explain me why the following program does not work?
Which pages of Unicode are implemented?
the u22xx and 2Axx are math symbols and extensions
Thanks
Heberto Ghezzo
Montreal

 plot(1)
 text(1.0,1.2,a \u2A8A b \u222C c \u5222 d, cex=2)
 text(1.0,1.1, \u222C , cex=2)
 text(1.0,1.0,b \u222C c, cex=2)
 text(1.2,1.0,c \u5222 d, cex=2)
 text(0.8,0.9,a \u2A8A b, cex=2)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Adding vertical line to histogram and qplot stacked plot

2009-01-29 Thread Jason Rupert

R-users it appears I am leaning on your knowledge once again.  Is there any way 
to add a vertical line to a histogram and qplot stacked plot?  Here is my 
current attempt:
 
qplot approach attempt:
qplot(Run, data = data_dataframe, breaks = breaks, fill = Temperature, main = 
short_title) + scale_x_continuous(Data) + scale_y_continuous(Freq) 
+  scale_fill_discrete(Temperature) + scale_fill_manual(values = c(LOW = 
blue, AMB =black, HIGH = red)) +  geom_abline(v = HighVal, col = 
dodgerblue3, lty=dotdash) 
 
 
hist approach attempt:
hist(data_dataframe, breaks = breaks,  col = dodgerblue3, xlab=Data, freq = 
TRUE, main=) +abline(v = HighVal, col = dodgerblue3, lty=dotdash)   
 
Neither seem to be working.  I think I am missing something small to get this 
to work.  
 
Thank you again for any feedback you can provide. 


  
[[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] Adding vertical line to histogram and qplot stacked plot

2009-01-29 Thread Felipe Carrillo

Jason:
Check Hadley's page, there's a few examples there. Good luck

http://had.co.nz/ggplot2/geom_vline.html

Felipe D. Carrillo  
Supervisory Fishery Biologist  
Department of the Interior  
US Fish  Wildlife Service  
California, USA


--- On Thu, 1/29/09, Jason Rupert jasonkrup...@yahoo.com wrote:

 From: Jason Rupert jasonkrup...@yahoo.com
 Subject: [R] Adding vertical line to histogram and qplot stacked plot
 To: R-help@r-project.org
 Date: Thursday, January 29, 2009, 11:03 AM
 R-users it appears I am leaning on your knowledge once
 again.  Is there any way to add a vertical line to a
 histogram and qplot stacked plot?  Here is my
 current attempt:
  
 qplot approach attempt:
 qplot(Run, data = data_dataframe, breaks = breaks, fill =
 Temperature, main = short_title)
 + scale_x_continuous(Data) +
 scale_y_continuous(Freq)
 +  scale_fill_discrete(Temperature) +
 scale_fill_manual(values = c(LOW = blue, AMB
 =black, HIGH = red)) + 
 geom_abline(v = HighVal, col = dodgerblue3,
 lty=dotdash) 
  
  
 hist approach attempt:
 hist(data_dataframe, breaks = breaks,  col =
 dodgerblue3, xlab=Data, freq = TRUE,
 main=) +abline(v = HighVal, col =
 dodgerblue3, lty=dotdash)   
  
 Neither seem to be working.  I think I am missing
 something small to get this to work.  
  
 Thank you again for any feedback you can provide. 
 
 
   
   [[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] Unicode

2009-01-29 Thread Duncan Murdoch


On 1/29/2009 1:57 PM, R Heberto Ghezzo, Dr wrote:

Hello, can somebody explain me why the following program does not work?


It works for me, though \u2A8A and \u222C display as boxes, something 
like the way they do in Firefox if I look at this page:


http://unicode.coeurlumiere.com/?n=8192

(but I don't see the embedded hex code.)


Which pages of Unicode are implemented?


Whatever your system and the graphic driver supports.  R isn't doing 
this by itself, it's asking the underlying system to do it.


Duncan Murdoch


the u22xx and 2Axx are math symbols and extensions
Thanks
Heberto Ghezzo
Montreal

 plot(1)
 text(1.0,1.2,a \u2A8A b \u222C c \u5222 d, cex=2)
 text(1.0,1.1, \u222C , cex=2)
 text(1.0,1.0,b \u222C c, cex=2)
 text(1.2,1.0,c \u5222 d, cex=2)
 text(0.8,0.9,a \u2A8A b, cex=2)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Adding vertical line to histogram and qplot stacked plot

2009-01-29 Thread Duncan Murdoch


On 1/29/2009 2:03 PM, Jason Rupert wrote:

R-users it appears I am leaning on your knowledge once again.  Is there any way to add a 
vertical line to a histogram and qplot stacked plot?  Here is my current 
attempt:
 
qplot approach attempt:
qplot(Run, data = data_dataframe, breaks = breaks, fill = Temperature, main = short_title) + scale_x_continuous(Data) + scale_y_continuous(Freq) +  scale_fill_discrete(Temperature) + scale_fill_manual(values = c(LOW = blue, AMB =black, HIGH = red)) +  geom_abline(v = HighVal, col = dodgerblue3, lty=dotdash) 
 
 
hist approach attempt:
hist(data_dataframe, breaks = breaks,  col = dodgerblue3, xlab=Data, freq = TRUE, main=) +abline(v = HighVal, col = dodgerblue3, lty=dotdash)   
 



Are you using ggplot2?  It uses that notation of adding elements to a 
plot.  Classic graphics doesn't, so you would just give the different 
commands separately on each line, e.g.


hist(data_dataframe, breaks = breaks,  col = dodgerblue3, xlab=Data, 
freq = TRUE, main=)


abline(v = HighVal, col = dodgerblue3, lty=dotdash)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] ggplot2 - how to change location / position of wind rose axis labels?

2009-01-29 Thread Darren Norris


Many thanks Hadley,
The solution you gave works well moving label locations within the area 
of the axis but does not quite provide what I was after. Apologies for 
the confusion - I've only been using ggplot2  it for a day.

Compare the bar chart (from previous example):
awind-doh + geom_bar(width = 1)

This shows x axis, gap beneath the x axis providing a separation to tick 
lines with x axis labels below tick lines.

The wind rose:
awind-doh + geom_bar(width = 1) + coord_polar()

Does not have the same separation between x axis and labels. How do I 
move all labels outside of the axis area (e.g. move label 8 to the left, 
3 to the right etc on the wind rose above)?

Thanks again for all your support,
Darren





hadley wickham wrote:
 Hi Darren,

 Sorry for taking so long to get back to you - my hard drive died and
 it's taking me a while to get up and running again of backups etc.

 One solution to your problem is to draw the axis labels yourself:

 labels - data.frame(rrating = 1:10, count = 19000)

 awind +
opts(axis.text.x = theme_blank()) +
geom_text(aes(label = rrating, y = count, fill = NULL), data = labels)

 Adjusting count as needed to get them in the right position.  I'll
 think about how to do this better in general.

 Regards,

 Hadley

 On Fri, Jan 23, 2009 at 4:39 PM, Darren Norris  wrote:
   
 Dear R users,
 First just want to say thank you to all for developing such a wonderful
 software and packages.

 I need to produce a wind rose plot. Tried with packages circular and plotrix
 and couldn't quite get what I want. Moved to package ggplot2 and it's going
 great. However stuck in how to move axis labels.

 I am using the wind rose from the help to learn how to do what I need (code
 example below). The plot produced has axis labels on top of the axis line
 (not sure if it's actually the plot or axis line...) - which means the
 labels are unclear for what I need to produce.

 How can I move the labels to be outside of the line? I have read the online
 book ( http://had.co.nz/ggplot2/book/ ), help, tried changing various theme
 and scale settings and searched the package website and mailing list (
 http://had.co.nz/ggplot2/ ). If the answer is there I'm too stupid to see it
 - do I need to play with grobs? If so how?

 Any help much appreciated,
 Darren

 R version 2.8.1 (2008-12-22)
 i386-pc-mingw32

 locale:
 LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
 Kingdom.1252;LC_MONETARY=English_United
 Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

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

 other attached packages:
 [1] ggplot2_0.8.1 reshape_0.8.2 plyr_0.1.4proto_0.3-8


 #Using an example form the coord_polar help
 library(ggplot2)
 movies$rrating - factor(round_any(movies$rating, 1))
 movies$budgetq - factor(chop(movies$budget, 4), labels = 1:4)
 doh - ggplot(movies, aes(x = rrating, fill = budgetq))
 doh + geom_bar(width = 1) + coord_polar()

 #Now with my theme (hacked from theme_bw) getting close to what I need

 theme_bwdn-function (base_size = 12)
 {
structure(list(axis.line = theme_blank(), axis.text.x = theme_text(size
 = base_size *
   1, lineheight = 0.9, vjust = 1), axis.text.y = theme_blank(),
 axis.ticks = theme_segment(colour = black,
size = 0.2), axis.title.x = theme_blank(), axis.title.y =
 theme_blank(), axis.ticks.length = unit(0,
lines), axis.ticks.margin = unit(0, lines), legend.background =
 theme_rect(colour = NA),
legend.key = theme_rect(colour = grey80), legend.key.size =
 unit(1.2,
lines), legend.text = theme_text(size = base_size *
0.8), legend.title = theme_text(size = base_size *
0.8, face = bold, hjust = 0), legend.position = right,
panel.background = theme_rect(fill = white, colour = NA),
panel.border = theme_rect(fill = NA, colour = grey50),
panel.grid.major = theme_line(colour = black, size = 0.2),
panel.grid.minor = theme_line(colour = black, size = 0.5),
panel.margin = unit(0.25, lines), strip.background =
 theme_rect(fill = grey80,
colour = grey50), strip.label = function(variable,
value) value, strip.text.x = theme_text(size = base_size *
0.8), strip.text.y = theme_text(size = base_size *
0.8, angle = -90), plot.background = theme_rect(colour = NA),
plot.title = theme_text(size = base_size * 1.2), plot.margin =
 unit(rep(0.5,
4), lines)), class = options)
 }

 awind-doh + geom_bar(width = 1) + coord_polar()

 #produces the wind rose but how to move axis labels?
 awind + theme_bwdn()




 --
 View this message in context: 
 http://www.nabble.com/ggplot2---how-to-change-location---position-of-wind-rose-axis-labels--tp21635526p21635526.html
 Sent from the R help mailing list archive at Nabble.com.

 __
 R-help@r-project.org mailing

[R] Importing data from clipboard on Mac OSX

2009-01-29 Thread Christian Anderson

Hello R-Help,

I noticed that there is a thread about importing data from the clipboard
that is very poorly answered in the forum. One user suggests giving up, the
other gives a solution that echoes the clipboard, but that's exactly the
same as just ctrl-p. As I am asked this question at least once a week in my
very small home institution, I'm sure many people want to know. If you could
post somewhere that for most Macs you can 

read.table(pipe(pbpaste)) will work for almost all applications, that
would be very helpful.

Thank you,
Christian Anderson

 


[[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] Welcome to the R-help mailing list

2009-01-29 Thread Bishwa S Koirala


On Thu, 29 Jan 2009 20:55:19 +0100
 r-help-requ...@r-project.org wrote:

Welcome to the R-help@r-project.org mailing list!

To post to this list, send your email to:

 r-help@r-project.org

General information about the mailing list is at:

 https://stat.ethz.ch/mailman/listinfo/r-help

If you ever want to unsubscribe or change your options 
(eg, switch to
or from digest mode, change your password, etc.), visit 
your

subscription page at:

 https://stat.ethz.ch/mailman/options/r-help/bkoirala%40unm.edu

You can also make such adjustments via email by sending 
a message to:


 r-help-requ...@r-project.org

with the word `help' in the subject or body (don't 
include the
quotes), and you will get back a message with 
instructions.


You must know your password to change your options 
(including changing

the password, itself) or to unsubscribe.  It is:

 bish90nu

Normally, Mailman will remind you of your r-project.org 
mailing list
passwords once every month, although you can disable 
this if you
prefer.  This reminder will also include instructions on 
how to
unsubscribe or change your account options.  There is 
also a button on
your options page that will email your current password 
to you.


Bishwa S Koirala
Department of Economics
University of New Mexico
Albuquerque, NM 87131

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] tab characters

2009-01-29 Thread Nick Matzke


Hi all,

Working at the R command line, how do I get strings to display e.g. tab 
or newline characters as they should be displayed, rather than as e.g. 
\n or \t?


e.g.:
 x=\t
 x=\t
 x
[1] \t
 print(x)
[1] \t

--

Nicholas J. Matzke
Ph.D. student, Graduate Student Researcher
Huelsenbeck Lab
Center for Theoretical Evolutionary Genomics
4151 VLSB (Valley Life Sciences Building)
Department of Integrative Biology
University of California, Berkeley

Lab websites:
http://ib.berkeley.edu/people/lab_detail.php?lab=54
http://fisher.berkeley.edu/cteg/hlab.html
Dept. personal page: 
http://ib.berkeley.edu/people/students/person_detail.php?person=370

Lab personal page: http://fisher.berkeley.edu/cteg/members/matzke.html
Lab phone: 510-643-6299
Dept. fax: 510-643-6264
Cell phone: 510-301-0179
Email: mat...@berkeley.edu

Mailing address:
Department of Integrative Biology
3060 VLSB #3140
Berkeley, CA 94720-3140

-
[W]hen people thought the earth was flat, they were wrong. When people 
thought the earth was spherical, they were wrong. But if you think that 
thinking the earth is spherical is just as wrong as thinking the earth 
is flat, then your view is wronger than both of them put together.


Isaac Asimov (1989). The Relativity of Wrong. The Skeptical Inquirer, 
14(1), 35-44. Fall 1989.

http://chem.tufts.edu/AnswersInScience/RelativityofWrong.htm

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

Re: [R] Stacking data

2009-01-29 Thread hadley wickham

 Â Â Â Â Â Â Â Â Â Â  Â Â  27193Â Â Â  30949 Â  Â Â  V2

 But i have 41 columns (age columnÂ  +Â  40 individuals)
 I have the following script but an error is thrown up
 can anyone help, where am i going wrong

 zz - read.csv(Filename.csv,strip.white = TRUE)

Try this:

install.packages(reshape)
library(reshape)
zzz - melt(zz)

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] RDieHarder with file input invocation

2009-01-29 Thread Tomek S.

I got hold of the 0.0.7 RDieHarder binary package as zip for Windows from 
CRANextras but don't know how to invoke its main method dieharder() in R  with 
file_input as RNG argument - how to specify path to that input file ? 
the 0.0.7 as opposed to 0.1.0 doesn't seem to have ''inputfile'' arg for path 
specification. So is file_input prior to 0.1.0 verison possible ? If so then 
how and does the input file have to be in a special format, have any header 
info etc to be processed properly by that dieharder() method. 
Thanks in advance.
Tom


[[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] assign point values as labels in spplot

2009-01-29 Thread emorway


In the code to follow, I'm trying to label points with their corresponding
values but have been unsuccessful after many variations to the following
code.  The code below creates the plot I want, I simply cannot get the black
points (+) to display the actual value.  I'm guessing the problem is
somewhere in the second to last line of code (starts with pts-...).  I
have attached the two text files needed to run the code.  
http://www.nabble.com/file/p21734824/R_Hk_Krig_File_Log.txt
R_Hk_Krig_File_Log.txt 
http://www.nabble.com/file/p21734824/Predict_Location_XY.txt
Predict_Location_XY.txt 
Eric

library(gstat)

K.dat-read.table(C:/temp/R_Hk_Krig_File_Log.txt, header = TRUE)
my.dat-data.frame(K.dat)
attach(my.dat)
coordinates(my.dat)=~X+Y

pred.dat-read.table(C:/temp/Predict_Location_XY.txt, header = FALSE)
names(pred.dat)-c(x,y,p)
my.pred.dat-data.frame(pred.dat)
coordinates(my.pred.dat)-~x+y
gridded(my.pred.dat) = TRUE

#Set up 2 exponential variograms
lzm.vgm.exp.noNug-vgm(1.1,Exp,2000,0)
lzm.vgm.exp.Nug-vgm(0.7,Exp,3000,0.4)

#Krige the exponential variograms
lzm.krig.exp.noNug-krige(Kh_Log10~1,my.dat,my.pred.dat,model=lzm.vgm.exp.noNug)
lzm.krig.exp.Nug-krige(Kh_Log10~1,my.dat,my.pred.dat,model=lzm.vgm.exp.Nug)

lzm.krig.exp.noNug$var1.pred-10^lzm.krig.exp.noNug$var1.pred
lzm.krig.exp.Nug$var1.pred-10^lzm.krig.exp.Nug$var1.pred

library(sp)
library(lattice)
pts-list(sp.points,K.dat,pch=3,col=black,labels=as.character(K.dat$Kh))
spplot(lzm.krig.exp.Nug[var1.pred], scales=list(draw=TRUE), 
xlab=Easting,ylab=Northing,cuts=25,key.space=right,cex=1.1,
col.regions=terrain.colors(25),main=Hydraulic Conductivity of Layer 2,
sp.layout=list(pts))

-- 
View this message in context: 
http://www.nabble.com/assign-point-values-as-labels-in-spplot-tp21734824p21734824.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] Welcome to the R-help mailing list

2009-01-29 Thread Johannes Huesing

Bishwa S Koirala bkoir...@unm.edu [Thu, Jan 29, 2009 at 08:56:17PM CET]:
 You must know your password to change your options (including changing
 the password, itself) or to unsubscribe.  It is:


It is definitely a good idea to change your password ASAP.


-- 
Johannes Hüsing   There is something fascinating about science. 
  One gets such wholesale returns of conjecture 
mailto:johan...@huesing.name  from such a trifling investment of fact.  
  
http://derwisch.wikidot.com (Mark Twain, Life on the Mississippi)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] tab characters

2009-01-29 Thread jim holtman

?cat

 x - '\t'
 print(x)
[1] \t
 cat(x)



On Thu, Jan 29, 2009 at 3:13 PM, Nick Matzke mat...@berkeley.edu wrote:
 Hi all,

 Working at the R command line, how do I get strings to display e.g. tab or
 newline characters as they should be displayed, rather than as e.g. \n or
 \t?

 e.g.:
 x=\t
 x=\t
 x
 [1] \t
 print(x)
 [1] \t

 --
 
 Nicholas J. Matzke
 Ph.D. student, Graduate Student Researcher
 Huelsenbeck Lab
 Center for Theoretical Evolutionary Genomics
 4151 VLSB (Valley Life Sciences Building)
 Department of Integrative Biology
 University of California, Berkeley

 Lab websites:
 http://ib.berkeley.edu/people/lab_detail.php?lab=54
 http://fisher.berkeley.edu/cteg/hlab.html
 Dept. personal page:
 http://ib.berkeley.edu/people/students/person_detail.php?person=370
 Lab personal page: http://fisher.berkeley.edu/cteg/members/matzke.html
 Lab phone: 510-643-6299
 Dept. fax: 510-643-6264
 Cell phone: 510-301-0179
 Email: mat...@berkeley.edu

 Mailing address:
 Department of Integrative Biology
 3060 VLSB #3140
 Berkeley, CA 94720-3140

 -
 [W]hen people thought the earth was flat, they were wrong. When people
 thought the earth was spherical, they were wrong. But if you think that
 thinking the earth is spherical is just as wrong as thinking the earth is
 flat, then your view is wronger than both of them put together.

 Isaac Asimov (1989). The Relativity of Wrong. The Skeptical Inquirer,
 14(1), 35-44. Fall 1989.
 http://chem.tufts.edu/AnswersInScience/RelativityofWrong.htm

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




-- 
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] Graphic device graphics primitives

2009-01-29 Thread Greg Snow

Sigbert,

The plot2script function in the TeachingDemos package does essentially what 
Duncan talks about for you.  Create your plot then run the function giving it a 
filename to save the info into (or run without arguments and then past into a 
script window or text editor (only tested on windows, if this does not work, go 
with the filename option)).  The same warnings (and possibly more) apply, but 
this lets you see the steps to recreate the plot (and you can try modifying 
some parts and running the script as a way to update the plot).

One thing this tends to show is the number and type of steps that the plotting 
process gets broken down into.  

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Duncan Murdoch
 Sent: Thursday, January 29, 2009 4:25 AM
 To: Sigbert Klinke
 Cc: r-help@r-project.org
 Subject: Re: [R] Graphic device  graphics primitives
 
 Sigbert Klinke wrote:
  Hi,
 
  I know that some graphics devices in R store graphics primitives such
  that a redraw is possible (e.g. when resizing the window). Is it
  possible to get the current number of stored graphic primitives?
 
  Thanks in advance
 
  Sigbert Klinke
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 You could examine the results of recordPlot, but note the warnings that
 the format is not guaranteed to be stable across R versions.
 
 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.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread SnowManPaddington


Hi ya, thanks a lot everyone!! I changed rr:ii-1 to rr:(ii-1) and the code
works!!! I finally get some estimates from the optimization function (i am
doing a logit model with 2 segments). Thanks thanks!!!
I didn't realize rr:(ii-1) and rr:ii-1 would make such a big difference,
especially because the professor used rr:ii-1 in his Gauss code. I didn't
realize it means so much different in R!



davidr-2 wrote:
 
 The lines below made me understand clearly. Maybe they are already in
 some documentation,
 but if not, it might help others to avoid my misunderstanding.
 Thanks to all for the clarifications.
 -- David
 
 
 -Original Message-
 From: Duncan Murdoch [mailto:murd...@stats.uwo.ca] 
 Sent: Thursday, January 29, 2009 9:54 AM
 To: David Reiner dav...@rhotrading.com
 Cc: Henrik Bengtsson; r-help@r-project.org; SnowManPaddington
 Subject: [SPAM] - Re: [R] scoping rules for 'for' [was - Re: for/if loop
 ] - Found word(s) list error in the Text body
 
 snip
   Feel free to do what you like to the variable within the 
 loop (though you might cause Luke to grind his teeth when it messes up 
 his compiler).  It will be set to the next value in the 1:10 vector the 
 next time it comes back to the top.
 
 snip
 
 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.
 
 

-- 
View this message in context: 
http://www.nabble.com/for-if-loop-tp21701496p21736578.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] bootstrapping in regression

2009-01-29 Thread Tom Backer Johnsen


Tom Backer Johnsen wrote:

Thomas Mang wrote:

Hi,

Please apologize if my questions sounds somewhat 'stupid' to the 
trained and experienced statisticians of you. Also I am not sure if I 
used all terms correctly, if not then corrections are welcome.


I have asked myself the following question regarding bootstrapping in 
regression:
Say for whatever reason one does not want to take the p-values for 
regression coefficients from the established test statistics 
distributions (t-distr for individual coefficients, F-values for 
whole-model-comparisons), but instead apply a more robust approach by 
bootstrapping.


In the simple linear regression case, one possibility is to randomly 
rearrange the X/Y data pairs, estimate the model and take the 
beta1-coefficient. Do this many many times, and so derive the null 
distribution for beta1. Finally compare beta1 for the observed data 
against this null-distribution.


There is a very basic difference between bootstrapping and random 
permutations.  What you are suggesting is to shuffle values between 
cases or rows in the frame.  That amounts to a variant of a permutation 
test, not a bootstrap.


What you do in a bootstrap test is different, you regard your sample as 
a population and then sample from that population (with replacement), 
normally by extracting a large number of random samples of the same size 
as the original sample and do the computations for whatever you are 
interested in for each sample.


In other words, with bootstrapping, the pattern of values within each 
case or row is unchanged, and you sample complete cases or rows.  With a 
permutation test you keep the original sample of cases or rows, but 
shuffle the observations on the same variable between cases or rows.


Have a look at the 'boot' package.

Tom


What I now wonder is how the situation looks like in the multiple 
regression case. Assume there are two predictors, X1 and X2. Is it 
then possible to do the same, but just only rearranging the values of 
one predictor (the one of interest) at a time? Say I want again to 
test beta1. Is it then valid to many times randomly rearrange the X1 
data (and keeping Y and X2 as observed), fit the model, take the beta1 
coefficient, and finally compare the beta1 of the observed data 
against the distributions of these beta1s ?
For X2, do the same, randomly rearrange X2 all the time while keeping 
Y and X1 as observed etc.

Is this valid ?

Second, if this is valid for the 'normal', fixed-effects only 
regression, is it also valid to derive null distributions for the 
regression coefficients of the fixed effects in a mixed model this 
way? Or does the quite different parameters estimation calculation 
forbid this approach (Forbid in the sense of bogus outcome) ?


Thanks, Thomas

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

and provide commented, minimal, self-contained, reproducible code.






--
++
| Tom Backer Johnsen, Psychometrics Unit,  Faculty of Psychology |
| University of Bergen, Christies gt. 12, N-5015 Bergen,  NORWAY |
| Tel : +47-5558-9185Fax : +47-5558-9879 |
| Email : bac...@psych.uib.noURL : http://www.galton.uib.no/ |
++

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Optim error: initial value in 'vmmin' is not finite

2009-01-29 Thread SnowManPaddington


Error in optim(method = BFGS, c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  : 
  initial value in 'vmmin' is not finite

I am running a logit model with latent class segments. I successfully got
estimates for 2 segments. However when I tried to increase the no. of
segments, I got this error message at the end. I checked my code again but
can't find anything wrong. Is this error message related to something
different from the code? Is that anything to do with the initial value? what
does vmmin mean exactly?

Thanks a lot!!


-- 
View this message in context: 
http://www.nabble.com/Optim-error%3A-initial-value-in-%27vmmin%27-is-not-finite-tp21736737p21736737.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] Importing data from clipboard on Mac OSX

2009-01-29 Thread Charles C. Berry


On Thu, 29 Jan 2009, Christian Anderson wrote:


Hello R-Help,

I noticed that there is a thread about importing data from the clipboard
that is very poorly answered in the forum. One user suggests giving up, the
other gives a solution that echoes the clipboard, but that's exactly the
same as just ctrl-p. As I am asked this question at least once a week in my
very small home institution, I'm sure many people want to know. If you could
post somewhere that for most Macs you can

read.table(pipe(pbpaste)) will work for almost all applications, that
would be very helpful.



?pipe says:


 Clipboard:

 ...

 MacOS X users can use 'pipe(pbpaste)' and 'pipe(pbcopy, w)'
 to read from and write to that system's clipboard.


HTH,

Chuck




Thank you,
Christian Anderson




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



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

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

Re: [R] Goodness of fit for gamma distributions

2009-01-29 Thread Albyn Jones


it is easy to make a qqplot for the gamma; suppose that the sample parameters
are 1.101 and 2.49, the data in x:

 plot(qgamma(ppoints(x),1.101,2.49),sort(x))

see also lattice:qqmath

albyn

Quoting Dan31415 d.m.mitch...@reading.ac.uk:



Ah yes, that does produce a nice plot. Can i just ask what exactly it is
showing. It seems to me to be a sort of Q-Q plot but with a different set of
axes. Is this correct, if so do the same interpretation rules apply for this
plot, i.e. departures from either end of the curve show poor fitting of the
extreme data.

thanks for your help Remko, its been very helpful.

Dann



Remko Duursma-2 wrote:


It sounds like you just want to graph it though. For gammas, it's nice
to graph the log of the density, because
the tail is so thin and long, so you don't see much otherwise:

mydata - rgamma(1, shape=1.1, rate=2.5)

# now suppose you fit a gamma distribution, and get these estimated
parameters:
shapeest - 1.101
rateest - 2.49

h - hist(mydata, breaks=50, plot=FALSE)
plot(h$mids, log(h$density))
curve(log(dgamma(x, shape=shapeest, rate=rateest)), add=TRUE)


#Remko


-
Remko Duursma
Post-Doctoral Fellow

Centre for Plant and Food Science
University of Western Sydney
Hawkesbury Campus
Richmond NSW 2753

Dept of Biological Science
Macquarie University
North Ryde NSW 2109
Australia

Mobile: +61 (0)422 096908



On Wed, Jan 28, 2009 at 1:13 AM, Dan31415 d.m.mitch...@reading.ac.uk
wrote:


Thanks for that Remko, but im slightly confused because isnt this testing
the
goodness of fit of 2 slightly different gamma distributions, not of how
well
a gamma distribution is representing the data.

e.g.

data.vec-as.vector(data)

(do some mle to find the parameters of a gamma distribution for data.vec)

xrarea-seq(-2,9,0.05)
yrarea-dgamma(xrarea,shape=7.9862,rate=2.6621)

so now yrarea is the gamma distribution and i want to compare it with
data.vec to see how well it fits.

regards,
Dann


Remko Duursma-2 wrote:


Hi Dann,

there is probably a better way to do this, but this works anyway:

# your data
gamdat - rgamma(1, shape=1, rate=0.5)

# comparison to gamma:
gamsam - rgamma(1, shape=1, rate=0.6)

qqplot(gamsam,gamdat)
abline(0,1)


greetings
Remko


-
Remko Duursma
Post-Doctoral Fellow

Centre for Plant and Food Science
University of Western Sydney
Hawkesbury Campus
Richmond NSW 2753

Dept of Biological Science
Macquarie University
North Ryde NSW 2109
Australia

Mobile: +61 (0)422 096908



On Tue, Jan 27, 2009 at 3:38 AM, Dan31415 d.m.mitch...@reading.ac.uk
wrote:


I'm looking for goodness of fit tests for gamma distributions with
large
data
sizes. I have a matrix with around 10,000 data values in it and i have
fitted a gamma distribution over a histogram of the data.

The problem is testing how well that distribution fits. Chi-squared
seems
to
be used more for discrete distributions and kolmogorov-smirnov seems
that
large sample sizes make it had to evaluate the D statistic. Also i
haven't
found a qq plot for gamma, although i think this might be an
appropriate
test.

in summary
-is there a gamma goodness of fit test that doesnt depend on the sample
size?
-is there a way of using qqplot for gamma distributions, if so how
would
you
calculate it from a matrix of data values?

regards,
Dann
--
View this message in context:
http://www.nabble.com/Goodness-of-fit-for-gamma-distributions-tp21668711p21668711.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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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://www.nabble.com/Goodness-of-fit-for-gamma-distributions-tp21668711p21686095.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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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://www.nabble.com/Goodness-of-fit-for-gamma-distributions-tp21668711p21725468.html

Sent from the R help mailing list

Re: [R] bootstrapping in regression

2009-01-29 Thread Greg Snow

What you are describing is actually a permutation test rather than a bootstrap 
(related concepts but with a subtle but important difference).

The way to do a permutation test with multiple x's is to fit the reduced model 
(use all x's other than x1 if you want to test x1) on the original data and 
store the fitted values and the residuals.

Permute the residuals (randomize their order) and add them back to the fitted 
values and fit the full model (including x1 this time) to the permuted data 
set.  Do this a bunch of times and it will give you the sampling distribution 
for the slope on x1 (or whatever your set of interest is) when the null 
hypothesis that it is 0 given the other variables in the model is true.

Permuting just x1 only works if x1 is orthogonal to all the other predictors, 
otherwise the permuting destroys the relationship with the other predictors and 
does not do the test you want.

Bootstrapping depends on sampling with replacement, not permuting, and is used 
more for confidence intervals than for tests (the reference by John Fox given 
to you in another reply can help if that is the approach you want to take).

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.s...@imail.org
801.408.8111


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Thomas Mang
 Sent: Thursday, January 29, 2009 9:44 AM
 To: r-h...@stat.math.ethz.ch
 Subject: [R] bootstrapping in regression
 
 Hi,
 
 Please apologize if my questions sounds somewhat 'stupid' to the
 trained
 and experienced statisticians of you. Also I am not sure if I used all
 terms correctly, if not then corrections are welcome.
 
 I have asked myself the following question regarding bootstrapping in
 regression:
 Say for whatever reason one does not want to take the p-values for
 regression coefficients from the established test statistics
 distributions (t-distr for individual coefficients, F-values for
 whole-model-comparisons), but instead apply a more robust approach by
 bootstrapping.
 
 In the simple linear regression case, one possibility is to randomly
 rearrange the X/Y data pairs, estimate the model and take the
 beta1-coefficient. Do this many many times, and so derive the null
 distribution for beta1. Finally compare beta1 for the observed data
 against this null-distribution.
 
 What I now wonder is how the situation looks like in the multiple
 regression case. Assume there are two predictors, X1 and X2. Is it then
 possible to do the same, but just only rearranging the values of one
 predictor (the one of interest) at a time? Say I want again to test
 beta1. Is it then valid to many times randomly rearrange the X1 data
 (and keeping Y and X2 as observed), fit the model, take the beta1
 coefficient, and finally compare the beta1 of the observed data against
 the distributions of these beta1s ?
 For X2, do the same, randomly rearrange X2 all the time while keeping Y
 and X1 as observed etc.
 Is this valid ?
 
 Second, if this is valid for the 'normal', fixed-effects only
 regression, is it also valid to derive null distributions for the
 regression coefficients of the fixed effects in a mixed model this way?
 Or does the quite different parameters estimation calculation forbid
 this approach (Forbid in the sense of bogus outcome) ?
 
 Thanks, Thomas
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-
 guide.html
 and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Power analysis for MANOVA?

2009-01-29 Thread Adam D. I. Kramer


Thanks for the response, Stephan.

Really, I am trying to say, My result is insignificant, my effect sizes are
tiny, you may want to consider the possibility that there really are no
meaningful differences. Computing post-hoc power makes a bit stronger of a
claim in this setting.

My real goal in this case was to put a single line on a poster that says,
Significance using our estimates would require N observations, which is
larger than the population. I am trying to solve for N.  N in this case is
a sort of effect size.  In this case, it is indeed a simple transformation
of Pillai's V and the p-value for the study, and I do not intend to suggest
that it is anything more than that.  However, I believe that the latter
effect size makes a much more compelling case, given that a lot of people
(such as yourself) don't have much experience with Pillai's V.

--Adam

On Wed, 28 Jan 2009, Stephan Kolassa wrote:


Hi Adam,

first: I really don't know much about MANOVA, so I sadly can't help you 
without learning about it an Pillai's V... which I would be glad to do, but I 
really don't have the time right now. Sorry!


Second: you seem to be doing a kind of post-hoc power analysis, my result 
isn't significant, perhaps that's due to low power? Let's look at the power 
of my experiment! My impression is that post-hoc power analysis and its 
interpretation is, shall we say, not entirely accepted within the statistical 
community, see:


Hoenig, J. M.,  Heisey, D. M. (2001, February). The abuse of power: The 
pervasive fallacy of power calculations for data analysis. The American 
Statistician, 55 (1), 1-6


And this:
http://staff.pubhealth.ku.dk/~bxc/SDC-courses/power.pdf

However, I am sure that lots of people can discuss this more competently than 
me...


Best wishes
Stephan


Adam D. I. Kramer schrieb:


On Mon, 26 Jan 2009, Stephan Kolassa wrote:


My (and, judging from previous traffic on R-help about power analyses,
also some other people's) preferred approach is to simply simulate an
effect size you would like to detect a couple of thousand times, run your
proposed analysis and look how often you get significance.  In your simple
case, this should be quite easy.


I actually don't have much experience running monte-carlo designs like
this...so while I'd certainly prefer a bootstrapping method like this one,
simulating the effect size given my constraints isn't something I've done
before.

The MANOVA procedure takes 5 dependent variables, and determines what
combination of the variables best discriminates the two levels of my
independent variable...then the discrimination rate is represented in the
statistic (Pillai's V=.00019), which is then tested (F[5,18653] = 0.71). 
So

coming up with a set of constraints that would produce V=.00019 given my
data set doesn't quite sound trivial...so I'll go for the par library
reference mentioned earlier before I try this.  That said, if anyone can
refer me to a tool that will help me out (or an instruction manual for 
RNG),

I'd also be much obliged.

Many thanks,
Adam




HTH,
Stephan


Adam D. I. Kramer schrieb:

Hello,

I have searched and failed for a program or script or method to
conduct a power analysis for a MANOVA. My interest is a fairly simple 
case

of 5 dependent variables and a single two-level categorical predictor
(though the categories aren't balanced).

If anybody happens to know of a script that will do this in R, I'd
love to know of it! Otherwise, I'll see about writing one myself.

What I currently see is this, from help.search(power):

stats::power.anova.test
Power calculations for balanced one-way
analysis of variance tests
stats::power.prop.test
Power calculations two sample test for
proportions
stats::power.t.test Power calculations for one and two sample t
tests

Any references on power in MANOVA would also be helpful, though of
course I will do my own lit search for them myself.

Cordially,
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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Optim error: initial value in 'vmmin' is not finite

2009-01-29 Thread Ben Bolker

SnowManPaddington wiwiana at gmail.com writes:

 
 
 Error in optim(method = BFGS, c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  : 
   initial value in 'vmmin' is not finite
 
 I am running a logit model with latent class segments. I successfully got
 estimates for 2 segments. However when I tried to increase the no. of
 segments, I got this error message at the end. I checked my code again but
 can't find anything wrong. Is this error message related to something
 different from the code? Is that anything to do with the initial value? what
 does vmmin mean exactly?
 

   As the error message says, there is a problem computing
the initial values for the model fitting procedure. These
initial values may be hidden within the code you are using.

  Please read the posting guide -- we need a lot more information
before we can help you.  If at all possible provide a reproducible
example, and tell us what packages/functions you are using (and
the output of sessionInfo() ).

  good luck
   Ben Bolker

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] convergence problem gamm / lme

2009-01-29 Thread geert aarts


problem solved: 

As is probably mostly the case, a convergence problem lies in the specification 
of the model or the data itself.

Some information: I was trying to model the spatial distribution of fish of a 
particular age. The raw observations consisted of the number of individuals of 
a particular length. We were interested in modelling the number of fish of age 
A per hour: Y*a/(t*s) where Y is the count, a is the probability that a fish of 
length l belongs to age A, t is the haul duration and s is the subsample factor 
(when not all fish caught are measured). So I decided to use Y as the response 
variable and log(t*s/a) as the model offset. In some cases a was 0, leading to 
a model offset of infinity. So I decided to replace those zeros (6 
observations) by a small value, but apparently not large enough. This must have 
caused the convergence problems in gamm.

After some reflection, I believe that one should often include the offset as 
model weights as well. E.g. assume you treat the log(fishing duration) as a 
model offset. If the fishing duration is very small compared to the rest or 
even 0, you havenât really made an observation, so it should receive a very 
small or zero model weight. In case of zero weight you might as well remove the 
observation from the analysis. At least that is how I see it.

So the solutions that worked:
1) Y*a as the response variable and log(t*s) as offset
2) Replacing the zero aâs by not such a small value
3) Removing the rows with zero aâs, using Y as the response variable and 
log(t*s/a) as the offset and weights. 

I believe that option 3 is most elegant.

Unfortunately it turned out that nobody could answer the question, because 
nobody knew the details of the data. 

Nevertheless, Simon thanks a lot for your replies!

Cheers Geert





Geert,

Sorry for slow reply... I don't see any obvious problems with what you've 
done, so I guess it's the usual problem that PQL just doesn't *have* to 
converge, and the bit of extra flexibility of using a smooth is too much for 
it in this case. If you send me the data offline I can dig a little bit more 
if you like (I'll only use the data for this purpose etc. etc.) 

You are right that PQL does the same thing for Poisson and quasi-poisson. I 
don't think there is an easy way to use the values for the reduced dataset 
fit in the full dataset fitting, unfortunately. 

Another option is to use `gam' to fit the random effects. It'll be a bit slow 
with 70+ random effects, as you have, and it's a bit more work to set up, but 
it should converge. See ?gam.models which has some examples showing how to do 
this.

best,
Simon
 

On Thursday 29 January 2009 08:20, geert aarts wrote:
 Simon, thanks for your reply and your suggestions.

 I fitted the following glmm's

 gamm3-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=l
ist(code_tripnr=~1),family=poisson))

 Which worked OK (see summary below)

 I also fitted a model using quasipoisson, but that didn't help. I actually
 also thought that glmmPQL and gamm estimate the dispersion parameter and
 hence assumes a quasipoisson distribution, even if you specify poisson. Is
 that correct?

 Finally I tried fitting a model to less data, and sometimes gamm managed to
 converge (see summary below). So would it be possible to use the parameter
 estimates from the model fitted to less data as starting values for the
 gamm fitted to the full data set? Or do you have any other suggestions?

 Thanks.
 Cheers Geert






 gamm3-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=l
ist(code_tripnr=~1),f

 amily=poisson))



 iteration
 1

 iteration
 2

 iteration
 3

detach(Disc_age)

 summary(gamm3)

 Linear
 mixed-effects model fit by maximum likelihood

  Data: NULL

   AIC BIC logLik

NA  NA
 NA



 Random
 effects:

  Formula: ~1 | code_tripnr

 (Intercept) Residual

 StdDev:
 0.001391914 231.9744



 Variance
 function:

  Structure: fixed weights

  Formula: ~invwt

 Fixed
 effects: count ~ offset(offsetter) + poly(lon, 3) * poly(lat, 3)

 Value
 Std.Error   DF t-value p-value

 (Intercept)-1.582 11.96 2024 -0.13232174  0.8947

 poly(lon,
 3)1  -4.048   1397.33 2024 -0.00289673  0.9977

 poly(lon,
 3)2 -22.013699.71 2024 -0.03145996  0.9749

 poly(lon,
 3)3  -8.538593.87 2024 -0.01437683  0.9885

 poly(lat,
 3)1-109.624666.05 2024 -0.16458856  0.8693

 poly(lat,
 3)2-104.179381.37 2024 -0.27316977  0.7848

 poly(lat,
 3)3 -10.661221.93 2024 -0.04803585  0.9617

 poly(lon,
 3)1:poly(lat, 3)1  4290.737  61369.98 2024
 0.06991589  0.9443

 poly(lon,
 3)2:poly(lat, 3)1  1853.559  36835.63 2024
 0.05031972  0.9599

 poly(lon,
 3)3:poly(lat, 3)1  -240.521  25771.80 2024 -0.00933272  0.9926

 poly(lon,
 3)1:poly(lat, 3)2  2540.147  41378.38 2024
 0.06138826

[R] How do I get my IT department to bless R?

2009-01-29 Thread Daniel Viar

I currently use R at work under the radar, but there's a chance I
could loose that access.  I'd like to get our company to feel
comfortable with open source and R in particular.  Does anyone have
any experience with their company's IT department and management that
they would be willing to share?  How does one get an all Microsoft
shop on board with allowing users to user R?  I know about the recent
NY Times article and recent news.  I'm afraid I may need some case
studies or examples of what other companies have done.

Any help would be greatly appreciated.

Thanks
Dan Viar
Chesapeake, VA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] bootstrapping in regression

2009-01-29 Thread Thomas Mang


Greg Snow wrote:

What you are describing is actually a permutation test rather than a bootstrap 
(related concepts but with a subtle but important difference).

The way to do a permutation test with multiple x's is to fit the reduced model 
(use all x's other than x1 if you want to test x1) on the original data and 
store the fitted values and the residuals.

Permute the residuals (randomize their order) and add them back to the fitted 
values and fit the full model (including x1 this time) to the permuted data 
set.  Do this a bunch of times and it will give you the sampling distribution 
for the slope on x1 (or whatever your set of interest is) when the null 
hypothesis that it is 0 given the other variables in the model is true.


Hi,

Thanks to you and Tom for the correction regarding bootstrapping vs 
permutation, and to Chuck for the cool link. Yes of course I described a 
permutation.


I have a question here: I am not sure if I understand your 'fit the full 
model ... to the permuted data set'. Am I correct to suppose that once 
the residuals of the reduced-model fit have been permuted and added back 
to the fitted values, the values obtained this way (fitted + permuted 
residuals) now constitute the new y-values to which the full model is 
fitted? Is that correct ?

Do you know if this procedure is also valid for a mixed-effects model ?

thanks a lot,
Thomas



Permuting just x1 only works if x1 is orthogonal to all the other predictors, 
otherwise the permuting destroys the relationship with the other predictors and 
does not do the test you want.

Bootstrapping depends on sampling with replacement, not permuting, and is used 
more for confidence intervals than for tests (the reference by John Fox given 
to you in another reply can help if that is the approach you want to take).

Hope this helps,



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] scoping rules for 'for' [was - Re: for/if loop ]

2009-01-29 Thread Gabor Grothendieck

On Thu, Jan 29, 2009 at 12:05 PM, Henrik Bengtsson h...@stat.berkeley.edu 
wrote:
 PS. About the double-letter index (e.g. ii vs. i); A few years ago
 someone suggested me to use this, because it is much easier to search
 for 'ii' in an editor compared with a single-letter 'i'.  So true.  I
 made the move immediately.


Its an interesting idea but it is ugly and surely one could just use
the capabilities of the editor for searching.  e.g. in vim /\i\ will find
i as a word but not as part of word.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Multidimensional scalling

2009-01-29 Thread Tomek Wlodarski

Dear R developers and users!

I have calculated metric MDS by cmdscale from matrix of distances
(dissimilarities).
I would like to ask you how can I estimate how well this new mapping
represents characteristic features of my data set?
Thank you for any suggestions.
Best,

tomek

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 do I get my IT department to bless R?

2009-01-29 Thread Erik Iverson

This is a very broad question, and the answer is going to depend on your 
particular situation, which we are not privy to.


I'll say two things.  First, you should try to figure out why they would not 
want you to run R, so you can address those reasons specifically.  Second, you 
might take a particular problem that you deal with, and specifically write out 
how R can make it easier, cheaper, more efficient to solve than the current 
solution.


Are there really still all-MS shops around?

Daniel Viar wrote:

I currently use R at work under the radar, but there's a chance I
could loose that access.  I'd like to get our company to feel
comfortable with open source and R in particular.  Does anyone have
any experience with their company's IT department and management that
they would be willing to share?  How does one get an all Microsoft
shop on board with allowing users to user R?  I know about the recent
NY Times article and recent news.  I'm afraid I may need some case
studies or examples of what other companies have done.

Any help would be greatly appreciated.

Thanks
Dan Viar
Chesapeake, VA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] gls prediction using the correlation structure in nlme

2009-01-29 Thread Dr Carbon

On Wed, Jan 28, 2009 at 9:51 AM, Dr Carbon drcar...@gmail.com wrote:
 How does one coerce predict.gls to incorporate the fitted correlation
 structure from the gls object into predictions? In the example below
 the AR(1) process with phi=0.545 is not used with predict.gls. Is
 there another function that does this? I'm going to want to fit a few
 dozen models varying in order from AR(1) to AR(3) and would like to
 look at the fits with the correlation structure included.

 Thanks in advance.

 -JC

 PS I am including the package maintainers on this post - does this
 constitute a maintainer-specific question in r-help etiquette?

 # example
 set.seed(123)
 x - arima.sim(list(order = c(1,0,0), ar = 0.7), n = 100)
 y -x + arima.sim(list(order = c(1,0,0), ar = 0.7), n = 100)
 x - c(x)
 y - c(y)
 lm1 - lm(y~x)
 ar(residuals(lm1)) # indicates an ar1 model
 cs1 - corARMA(p=1)
 fm1 - gls(y~x,corr=cs1)
 summary(fm1)
 # get fits
 fits - predict(fm1)
 # use coef to get fits
 fits2 - coef(fm1)[1] + coef(fm1)[2] * x
 plot(fits,fits2)


I think this is the way to do this?
b0 -  coef(fm1)[1]
b1 - coef(fm1)[2]
p1 - intervals(fm1)$corStruct[2]
y[i] = b0 + p1*y[i-1] + b1*x[i]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Q about how to use Anova.mlm

2009-01-29 Thread pgseye


Hi,

Am newish to stats and R, so I certainly appreciate any help. Basically I
have 50 inidividuals whom I have 6 photos each of their optic nerve head. I
want to check that the orientation of the nerve head is consistent, ie the 6
replicates show minimal or preferably no rotation differences. I'll draw an
arbitrary line between some blood vessels (same reference in each set of
replicates) and determine an angle of deviation from the vertical and that
angle will be my dependent variable.

Subject   Replicate   Angle of Deviation
1  1   x
1  2   x
1  3   x
1  4   x
1  5   x
1  6   x
2  1   x
2  2   x
2  3   x
2  4   x
2  5   x
2  6   x
etc

I'm wanting to test for Sphericity (because I've read that you should - is
this routine in a repeated measures ANOVA?) and can see that Anova.mlm in
the CAR package offers this in addition to the alternative Greenhouse and
Feldt tests.

I just don't really know how to perform the test - can someone give me some
help.

Thank you.

Paul

Dept of Ophthalmology
Uni Melbourne, Australia
-- 
View this message in context: 
http://www.nabble.com/Q-about-how-to-use-Anova.mlm-tp21739443p21739443.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] Regression

2009-01-29 Thread Edwin Wibisono

Hello, I have problem
I have a, and b in regression
then I can't plot x,y with (and) a, b lines

can you help me ?

thx

[[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] Matlab inv() and R solve() differences

2009-01-29 Thread Joseph P Gray

I submit the following matrix to both MATLAB and R

x=  0.133 0.254 -0.214 0.116
0.254 0.623 -0.674 0.139
   -0.214 -0.674 0.910 0.011
0.116 0.139  0.011 0.180

MATLAB's inv(x) provides the following
 137.21 -50.68 -4.70 -46.42
-120.71  27.28 -8.94 62.19
-58.15   6.93  -7.89  36.94
  8.35   11.17 10.42 -14.82

R's solve(x) provides:
261.94 116.22 150.92 -267.78
116.22 344.30 286.68 -358.30
150.92 286.68 252.96 -334.09
-267.78 =358.30 -334.09 475.22

inv(x)*x = I(4)
and solve(x)%*%x = I(4)

Is there a way to obtain the MATLAB result in R?

Thanks for any help.


Pat Gray

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Regression

2009-01-29 Thread stephen sefick

Provide reproducible code.  And you will get an answer, probably.

On Thu, Jan 29, 2009 at 8:29 PM, Edwin Wibisono edwin...@gmail.com wrote:
 Hello, I have problem
 I have a, and b in regression
 then I can't plot x,y with (and) a, b lines

 can you help me ?

 thx

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




-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 do I get my IT department to bless R?

2009-01-29 Thread Jim Porzak

Yes, Erik, there are all MS shops around! Ours happens to be one.

However, I have absolutely no push back from IT on my use of R to do
marketing analytics. The trick, Dan, is to deliver relevant and
actionable results to the business. Your champions will stick up for
you when, and if, you get any push back from IT.

HTH,
Jim Porzak
TGN.com
San Francisco, CA
http://www.linkedin.com/in/jimporzak
use R! Group SF: http://ia.meetup.com/67/



On Thu, Jan 29, 2009 at 5:13 PM, Erik Iverson iver...@biostat.wisc.edu wrote:
 This is a very broad question, and the answer is going to depend on your
 particular situation, which we are not privy to.

 I'll say two things.  First, you should try to figure out why they would not
 want you to run R, so you can address those reasons specifically.  Second,
 you might take a particular problem that you deal with, and specifically
 write out how R can make it easier, cheaper, more efficient to solve than
 the current solution.

 Are there really still all-MS shops around?

 Daniel Viar wrote:

 I currently use R at work under the radar, but there's a chance I
 could loose that access.  I'd like to get our company to feel
 comfortable with open source and R in particular.  Does anyone have
 any experience with their company's IT department and management that
 they would be willing to share?  How does one get an all Microsoft
 shop on board with allowing users to user R?  I know about the recent
 NY Times article and recent news.  I'm afraid I may need some case
 studies or examples of what other companies have done.

 Any help would be greatly appreciated.

 Thanks
 Dan Viar
 Chesapeake, VA

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

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


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

Re: [R] Multidimensional scalling

2009-01-29 Thread stephen sefick

It depends on what you mean?  If you would like a goodness of fit of
your ordination to your distance matrix then this is doable and I
would suggest that you look at the labdsv  tutorial -
http://ecology.msu.montana.edu/labdsv/R/


On Thu, Jan 29, 2009 at 6:50 PM, Tomek Wlodarski
tomek.wlodar...@gmail.com wrote:
 Dear R developers and users!

 I have calculated metric MDS by cmdscale from matrix of distances
 (dissimilarities).
 I would like to ask you how can I estimate how well this new mapping
 represents characteristic features of my data set?
 Thank you for any suggestions.
 Best,

 tomek

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




-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

-K. Mullis

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Q about how to use Anova.mlm

2009-01-29 Thread John Fox

Dear Paul,

First, to fit a multivariate linear model to your data, you'll have to
rearrange the data from long format (with one observation per replicate)
to wide format (with one observation per subject). If your data are in the
data frame Data, then you'd do something like:

Wide - reshape(Data, v.names=Angle, idvar=Subject, timevar=Replicate,
direction=wide)

Then, with the data in wide format, fit a multivariate linear model with
just a constant:

mod - lm(cbind(Angle.1, Angle.2, Angle.3, Angle.4, Angle.5, Angle.6) ~ 1,
data=DavisThin)

Finally, use Anova() to get the tests:

idata - data.frame(Replicate=c(Angle.1, Angle.2, Angle.3, Angle.4,
Angle.5, Angle.6))
summary(Anova(mod, idata=idata, idesign=~Replicate))

If I understand correctly what you want, this should give it to you.

As well, since your design has no between-subject factors and only a single
within-subject factor, you could also use anova() [i.e., anova.mlm()] to get
the same results.

I hope this helps,
 John

--
John Fox, Professor
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 pgseye
 Sent: January-29-09 7:40 PM
 To: r-help@r-project.org
 Subject: [R] Q about how to use Anova.mlm
 
 
 Hi,
 
 Am newish to stats and R, so I certainly appreciate any help. Basically I
 have 50 inidividuals whom I have 6 photos each of their optic nerve head.
I
 want to check that the orientation of the nerve head is consistent, ie the
6
 replicates show minimal or preferably no rotation differences. I'll draw
an
 arbitrary line between some blood vessels (same reference in each set of
 replicates) and determine an angle of deviation from the vertical and that
 angle will be my dependent variable.
 
 Subject Replicate   Angle of Deviation
 11   x
 12   x
 13   x
 14   x
 15   x
 16   x
 21   x
 22   x
 23   x
 24   x
 25   x
 26   x
 etc
 
 I'm wanting to test for Sphericity (because I've read that you should - is
 this routine in a repeated measures ANOVA?) and can see that Anova.mlm in
 the CAR package offers this in addition to the alternative Greenhouse and
 Feldt tests.
 
 I just don't really know how to perform the test - can someone give me
some
 help.
 
 Thank you.
 
 Paul
 
 Dept of Ophthalmology
 Uni Melbourne, Australia
 --
 View this message in context: http://www.nabble.com/Q-about-how-to-use-
 Anova.mlm-tp21739443p21739443.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-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Matlab inv() and R solve() differences

2009-01-29 Thread Kingsford Jones

I suppose the solution is unstable because x is ill-conditioned:

 x
   [,1]   [,2]   [,3]  [,4]
[1,]  0.133  0.254 -0.214 0.116
[2,]  0.254  0.623 -0.674 0.139
[3,] -0.214 -0.674  0.910 0.011
[4,]  0.116  0.139  0.011 0.180
 cor(x)
   [,1]   [,2]   [,3]   [,4]
[1,]  1.000  0.9963557 -0.9883690  0.8548065
[2,]  0.9963557  1.000 -0.9976663  0.8084090
[3,] -0.9883690 -0.9976663  1.000 -0.7663847
[4,]  0.8548065  0.8084090 -0.7663847  1.000

 kappa(x)
[1] 2813.326

hth,

Kingsford Jones

On Thu, Jan 29, 2009 at 7:00 PM, Joseph P Gray jpg...@uwm.edu wrote:
 I submit the following matrix to both MATLAB and R

 x=  0.133 0.254 -0.214 0.116
0.254 0.623 -0.674 0.139
   -0.214 -0.674 0.910 0.011
0.116 0.139  0.011 0.180

 MATLAB's inv(x) provides the following
  137.21 -50.68 -4.70 -46.42
 -120.71  27.28 -8.94 62.19
 -58.15   6.93  -7.89  36.94
  8.35   11.17 10.42 -14.82

 R's solve(x) provides:
 261.94 116.22 150.92 -267.78
 116.22 344.30 286.68 -358.30
 150.92 286.68 252.96 -334.09
 -267.78 =358.30 -334.09 475.22

 inv(x)*x = I(4)
 and solve(x)%*%x = I(4)

 Is there a way to obtain the MATLAB result in R?

 Thanks for any help.


 Pat Gray

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] Q about how to use Anova.mlm

2009-01-29 Thread John Fox

Dear Paul,

I noticed a typo in my response and some poor formatting in the email
message; please see below:

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On
 Behalf Of John Fox
 Sent: January-29-09 9:23 PM
 To: 'pgseye'
 Cc: r-help@r-project.org
 Subject: Re: [R] Q about how to use Anova.mlm
 
 Dear Paul,
 
 First, to fit a multivariate linear model to your data, you'll have to
 rearrange the data from long format (with one observation per replicate)
 to wide format (with one observation per subject). If your data are in
the
 data frame Data, then you'd do something like:
 
 Wide - reshape(Data, v.names=Angle, idvar=Subject,
timevar=Replicate,
 direction=wide)
 
 Then, with the data in wide format, fit a multivariate linear model with
 just a constant:
 
 mod - lm(cbind(Angle.1, Angle.2, Angle.3, Angle.4, Angle.5, Angle.6) ~ 1,
 data=DavisThin)

The data argument should be data=Wide (I adapted the code from an example I
already had and neglected to change the argument).

 
 Finally, use Anova() to get the tests:
 
 idata - data.frame(Replicate=c(Angle.1, Angle.2, Angle.3,
Angle.4,
 Angle.5, Angle.6))
 summary(Anova(mod, idata=idata, idesign=~Replicate))

Note that this is indeed two separate commands; they were apparently run
together by my mailer.

Regards,
 John

 
 If I understand correctly what you want, this should give it to you.
 
 As well, since your design has no between-subject factors and only a
single
 within-subject factor, you could also use anova() [i.e., anova.mlm()] to
get
 the same results.
 
 I hope this helps,
  John
 
 --
 John Fox, Professor
 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 pgseye
  Sent: January-29-09 7:40 PM
  To: r-help@r-project.org
  Subject: [R] Q about how to use Anova.mlm
 
 
  Hi,
 
  Am newish to stats and R, so I certainly appreciate any help. Basically
I
  have 50 inidividuals whom I have 6 photos each of their optic nerve
head.
 I
  want to check that the orientation of the nerve head is consistent, ie
the
 6
  replicates show minimal or preferably no rotation differences. I'll draw
 an
  arbitrary line between some blood vessels (same reference in each set of
  replicates) and determine an angle of deviation from the vertical and
that
  angle will be my dependent variable.
 
  Subject   Replicate   Angle of Deviation
  1  1   x
  1  2   x
  1  3   x
  1  4   x
  1  5   x
  1  6   x
  2  1   x
  2  2   x
  2  3   x
  2  4   x
  2  5   x
  2  6   x
  etc
 
  I'm wanting to test for Sphericity (because I've read that you should -
is
  this routine in a repeated measures ANOVA?) and can see that Anova.mlm
in
  the CAR package offers this in addition to the alternative Greenhouse
and
  Feldt tests.
 
  I just don't really know how to perform the test - can someone give me
 some
  help.
 
  Thank you.
 
  Paul
 
  Dept of Ophthalmology
  Uni Melbourne, Australia
  --
  View this message in context: http://www.nabble.com/Q-about-how-to-use-
  Anova.mlm-tp21739443p21739443.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-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] How do I get my IT department to bless R?

2009-01-29 Thread Warren Young


Daniel Viar wrote:

I'd like to get our company to feel
comfortable with open source


Anyone still denying, here in 2009, that open source offers serious 
business value is a dinosaur, doomed to extinction.  Their cerebella 
have calcified.  The balance tipped a decade ago.


Just like the real dinosaurs, extinction will only be fast on a 
geological time scale.  Don't expect your job to evaporate next year 
because they won't use open source.  Just expect that over the coming 
decades to be routinely outcompeted by the mammals.


Chances are, your company actually has embraced open source in some way. 
 One facile argument is to ask if they use Google.  Yes?  Google uses 
Linux, MySQL, and yes, even R, so your company does too, if indirectly. 
 Likely, some bit of open source has crept into your actual operation 
elsewhere besides your little R enclave.



How does one get an all Microsoft
shop on board with allowing users to user R?


Proceed the same way you already are.

It is as Gandhi said: First they ignore you, then they laugh at you, 
then they fight you, then you win.


Every revolution in corporate IT happened this way, including 
Microsoft's own rise to dominance.  (Remember Big Blue?)


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 do I get my IT department to bless R?

2009-01-29 Thread David M Smith

On Thu, Jan 29, 2009 at 2:29 PM, Daniel Viar dan.v...@gmail.com wrote:
 I currently use R at work under the radar, but there's a chance I
 could loose that access.  I'd like to get our company to feel
 comfortable with open source and R in particular.  Does anyone have
 any experience with their company's IT department and management that
 they would be willing to share?  How does one get an all Microsoft
 shop on board with allowing users to user R?  I know about the recent
 NY Times article and recent news.  I'm afraid I may need some case
 studies or examples of what other companies have done.

In many cases your IT department will feel secure with R if there's a
company behind it to offer technical support and offer a throat to
choke if problems arise.  (Whether it's the former or the latter that
is more significant depends on the company.)  There are some companies
out there that offer support subscriptions to R, including the one I
work for.

If you work in a regulated environment (such as clinical pharma with
21CFR11, or finance with Sarbanes-Oxley), there may also be some
nervousness about whether R can be compliant.  It almost certainly is,
but it often needs to be validated in your own environment. I wrote
about this recently (from the perspective of FDA validation) at
http://blog.revolution-computing.com/2009/01/analyzing-clinical-trial-data-with-r.html
.

In many companies IT departments are getting comfortable with relying
on FOSS applications, but the real successes (Linux, Apache, MySQL,
...) have come when there's a commercial company to back up the open
source community.

# David Smith

--
David M Smith da...@revolution-computing.com
Director of Community, REvolution Computing www.revolution-computing.com
Tel: +1 (206) 577-4778 x3203 (Seattle, USA)

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Matlab inv() and R solve() differences

2009-01-29 Thread Berwin A Turlach

G'day all,

On Thu, 29 Jan 2009 19:24:40 -0700
Kingsford Jones kingsfordjo...@gmail.com wrote:

 I suppose the solution is unstable because x is ill-conditioned:

While, as you show, x is ill-conditioned, I do not believe that this is
serious enough to explain the differences that Pat sees between MATLAB
and R.  

In fact, on one of our lab machines MATLAB Version 7.5.0.342 (R2007b),
August 15, 2007, yields the following result:

 x=[0.133 0.254 -0.214 0.116; 0.254 0.623 -0.674 0.139; -0.214 -0.674  0.910 
 0.011 ; 0.116 0.139 0.011 0.180]

x =
0.13300.2540   -0.21400.1160
0.25400.6230   -0.67400.1390
   -0.2140   -0.67400.91000.0110
0.11600.13900.01100.1800

 inv(x)  

ans =
  261.9426  116.2219  150.9174 -267.7793
  116.2219  344.3029  286.6735 -358.2959
  150.9174  286.6735  252.9553 -334.0920
 -267.7793 -358.2959 -334.0920  475.2252

which is consistent with the output of R's solve().  

But the matrix x is ill-conditioned enough for small changes in the
precision of the entries leading to big differences in the calculated
inverse matrix.

My guess is that Pat was not completely truthful in the description of
what he did.  Presumably, x was not submitted to MATLAB in the given
form but calculated (to a higher precision) in MATLAB before being fed
to MATLAB's inv() command.  Then x was transferred with a lower
precision (less significant digits) to R and submitted there to R's
solve() function. 

Thus coming back to Pat's original question:

 On Thu, Jan 29, 2009 at 7:00 PM, Joseph P Gray jpg...@uwm.edu wrote:
  [...]
  Is there a way to obtain the MATLAB result in R?

Presumably yes.  Feed x in the same precision to R's solve() function
as you feed it to MATLAB's inv() function.

HTH.

Cheers,

Berwin

=== Full address =
Berwin A TurlachTel.: +65 6516 4416 (secr)
Dept of Statistics and Applied Probability+65 6516 6650 (self)
Faculty of Science  FAX : +65 6872 3919   
National University of Singapore 
6 Science Drive 2, Blk S16, Level 7  e-mail: sta...@nus.edu.sg
Singapore 117546http://www.stat.nus.edu.sg/~statba

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

1 2 >

1 - 100 of 108 matches

Mail list logo