Re: [R] How to calculate the spatial correlation of several files?

2012-12-04 Thread Jonsson
So where is the final correlation map
Can we write it:
to.write =
file(paste(C:\\Users\\aalyaari\\desktop\\corr1.bin,sep=),wb)

writeBin(as.double(results[[.f]]), to.write, size = 4)




--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-calculate-the-spatial-correlation-of-several-files-tp4651888p4652004.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] xlsx file read in R

2012-12-04 Thread Martin Studer
Hi Nico,

please let me know the details of sessionInfo() in R. Also, what version of
Java are you running? If you could post the output of java -version from
the command line that would be great. Note that XLConnect requires at least
Java 1.6.

Best regards,
Martin




--
View this message in context: 
http://r.789695.n4.nabble.com/xlsx-file-read-in-R-tp4651829p4652003.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] question about sum of (column) elements in R

2012-12-04 Thread T Bal
Hi,
I have the following data:

0   12
1   10
1   4
1   6
1   7
1   13
2   21
2   23
2   20
3   18
3   17
3   16
3   27
3   33
4   11
4   8
4   19
4   16
4   9


In this data file I would like to sum the numbers of second column which
belong to the same number in the first column.
So the output would be:

0 12
1 40
2 64
3 111
etc.

Thank you.

Kind regards,
T. Bal

[[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] question about sum of (column) elements in R

2012-12-04 Thread Pascal Oettli

Hi,

You can follow this example:

test - structure(list(V1 = c(0L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L), V2 = c(12L, 10L, 4L, 6L,
7L, 13L, 21L, 23L, 20L, 18L, 17L, 16L, 27L, 33L, 11L, 8L, 19L,
16L, 9L)), .Names = c(V1, V2), class = data.frame, row.names = 
c(NA,-19L))


tapply(grid1$V2, grid1$V1, sum)

  0   1   2   3   4
 12  40  64 111  63

HTH
Pascal


Le 04/12/2012 16:59, T Bal a écrit :

Hi,
I have the following data:

0   12
1   10
1   4
1   6
1   7
1   13
2   21
2   23
2   20
3   18
3   17
3   16
3   27
3   33
4   11
4   8
4   19
4   16
4   9


In this data file I would like to sum the numbers of second column which
belong to the same number in the first column.
So the output would be:

0 12
1 40
2 64
 3 111
 etc.

Thank you.

Kind regards,
T. Bal

[[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] problem with factor levels

2012-12-04 Thread Jeremy.Shearman
Hi
  I have a data.frame with 371,718 obs. of 12 variables (see below for
an str). My problem is with V1, a Factor w/ 93144 levels, there should
actually be 93994 levels. Each entry looks like:
comp[number]_c[number]_seq[number]
for example
comp215489_c0_seq40
R is grouping as though the last number is a decimal for some reason, in
other words comp215489_c0_seq40 and comp215489_c0_seq4 are considered to be
the same. My problem is that they are not the same so when I group by this
factor I am losing 800 levels.

Here is an str

'data.frame':   371718 obs. of  12 variables:
 $ V1 : Factor w/ 93144 levels comp10_c0_seq1,..: 92271 91685 29 30
1564 1564 1623 91700 91701 91848 ...
 $ V2 : Factor w/ 17162 levels gi|345842331|ref|NM_001244016.1|,..: 10119
10779 13210 13210 11522 8115 13079 14493 14493 15858 ...
 $ V3 : num  95.5 90.2 98.7 99.2 81.4 ...
 $ V4 : int  335 153 237 122 258 127 306 258 120 177 ...
 $ V5 : int  15 15 3 1 38 19 20 23 5 9 ...
 $ V6 : int  0 0 0 0 4 2 0 0 0 0 ...
 $ V7 : int  1 45 1 43 1 129 1 54 1 70 ...
 $ V8 : int  335 197 237 164 254 254 306 311 120 246 ...
 $ V9 : int  6866 18 3172 3438 67 122 3927 42 346 195 ...
 $ V10: int  7200 170 3408 3559 318 247 4232 299 465 19 ...
 $ V11: num  7e-155 2e-46 4e-125 2e-61 3e-24 ...
 $ V12: num  545 184 446 234 111 69.9 448 329 198 280 ..



--
View this message in context: 
http://r.789695.n4.nabble.com/problem-with-factor-levels-tp4652006.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] problem with factor levels

2012-12-04 Thread Milan Bouchet-Valat
Le mardi 04 décembre 2012 à 00:34 -0800, Jeremy.Shearman a écrit :
 Hi
   I have a data.frame with 371,718 obs. of 12 variables (see below for
 an str). My problem is with V1, a Factor w/ 93144 levels, there should
 actually be 93994 levels. Each entry looks like:
 comp[number]_c[number]_seq[number]
 for example
 comp215489_c0_seq40
 R is grouping as though the last number is a decimal for some reason, in
 other words comp215489_c0_seq40 and comp215489_c0_seq4 are considered to be
 the same. My problem is that they are not the same so when I group by this
 factor I am losing 800 levels.
What format is your original data using? How do you import it?

Please provide us with an excerpt of your original file showing at least
two different values of V1 that are considered the same once imported in
R (which sounds very unlikely to me...).


Regards

 Here is an str
 
 'data.frame': 371718 obs. of  12 variables:
  $ V1 : Factor w/ 93144 levels comp10_c0_seq1,..: 92271 91685 29 30
 1564 1564 1623 91700 91701 91848 ...
  $ V2 : Factor w/ 17162 levels gi|345842331|ref|NM_001244016.1|,..: 10119
 10779 13210 13210 11522 8115 13079 14493 14493 15858 ...
  $ V3 : num  95.5 90.2 98.7 99.2 81.4 ...
  $ V4 : int  335 153 237 122 258 127 306 258 120 177 ...
  $ V5 : int  15 15 3 1 38 19 20 23 5 9 ...
  $ V6 : int  0 0 0 0 4 2 0 0 0 0 ...
  $ V7 : int  1 45 1 43 1 129 1 54 1 70 ...
  $ V8 : int  335 197 237 164 254 254 306 311 120 246 ...
  $ V9 : int  6866 18 3172 3438 67 122 3927 42 346 195 ...
  $ V10: int  7200 170 3408 3559 318 247 4232 299 465 19 ...
  $ V11: num  7e-155 2e-46 4e-125 2e-61 3e-24 ...
  $ V12: num  545 184 446 234 111 69.9 448 329 198 280 ..
 
 
 
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/problem-with-factor-levels-tp4652006.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] question about sum of (column) elements in R

2012-12-04 Thread Gerrit Eichner

Hi, T. Bal,

homework? Take a look at

?tapply


 Regards  --  Gerrit


On Tue, 4 Dec 2012, T Bal wrote:


Hi,
I have the following data:

0   12
1   10
1   4
1   6
1   7
1   13
2   21
2   23
2   20
3   18
3   17
3   16
3   27
3   33
4   11
4   8
4   19
4   16
4   9


In this data file I would like to sum the numbers of second column which
belong to the same number in the first column.
So the output would be:

0 12
1 40
2 64
   3 111
   etc.

Thank you.

Kind regards,
T. Bal

[[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] partial analisys of a time series

2012-12-04 Thread Antonio Silva
Dear list members

I want to analyze separately the months of a time series. In other words, I
want to plot and fit models for each month separately.

Taking the example of
http://a-little-book-of-r-for-time-series.readthedocs.org/en/latest/src/timeseries.html

births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
birthstimeseries - ts(births, frequency=12, start=c(1946,1))
birthstimeseries
plot.ts(birthstimeseries)
birthstimeseriesHW - HoltWinters(birthstimeseries)
plot(birthstimeseriesHW)

How to proceed the plotting and HoltWinters smoothing considereing only
Januarys, Februarys, etc. separately.

Thanks in advance.

Antonio Olinto

[[alternative HTML version deleted]]

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


Re: [R] problem with factor levels

2012-12-04 Thread PIKAL Petr
Hi

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Jeremy.Shearman
 Sent: Tuesday, December 04, 2012 9:35 AM
 To: r-help@r-project.org
 Subject: [R] problem with factor levels
 
 Hi
   I have a data.frame with 371,718 obs. of 12 variables (see below
 for an str). My problem is with V1, a Factor w/ 93144 levels, there
 should actually be 93994 levels. Each entry looks like:
 comp[number]_c[number]_seq[number]
 for example
 comp215489_c0_seq40
 R is grouping as though the last number is a decimal for some reason,
 in other words comp215489_c0_seq40 and comp215489_c0_seq4 are
 considered to be the same. My problem is that they are not the same so
 when I group by this factor I am losing 800 levels.
 

Hm. How did you constructed those factors?

 factor(c(comp215489_c0_seq40, comp215489_c0_seq4) )
[1] comp215489_c0_seq40 comp215489_c0_seq4 
Levels: comp215489_c0_seq4 comp215489_c0_seq40

gives me 2 levels as expected. I also doubt that R will do such stripping 
during reading from other file.

Regards
Petr


 Here is an str
 
 'data.frame': 371718 obs. of  12 variables:
  $ V1 : Factor w/ 93144 levels comp10_c0_seq1,..: 92271 91685 29
 30
 1564 1564 1623 91700 91701 91848 ...
  $ V2 : Factor w/ 17162 levels gi|345842331|ref|NM_001244016.1|,..:
 10119
 10779 13210 13210 11522 8115 13079 14493 14493 15858 ...
  $ V3 : num  95.5 90.2 98.7 99.2 81.4 ...
  $ V4 : int  335 153 237 122 258 127 306 258 120 177 ...
  $ V5 : int  15 15 3 1 38 19 20 23 5 9 ...
  $ V6 : int  0 0 0 0 4 2 0 0 0 0 ...
  $ V7 : int  1 45 1 43 1 129 1 54 1 70 ...
  $ V8 : int  335 197 237 164 254 254 306 311 120 246 ...
  $ V9 : int  6866 18 3172 3438 67 122 3927 42 346 195 ...
  $ V10: int  7200 170 3408 3559 318 247 4232 299 465 19 ...
  $ V11: num  7e-155 2e-46 4e-125 2e-61 3e-24 ...
  $ V12: num  545 184 446 234 111 69.9 448 329 198 280 ..
 
 
 
 --
 View this message in context: http://r.789695.n4.nabble.com/problem-
 with-factor-levels-tp4652006.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] problem with factor levels

2012-12-04 Thread Jeremy.Shearman
Oh, your skepticism was spot on!
I was using excel to check the output (silly, but I am still in the process
of moving from excel to R) and there was a discrepancy in the number of
output from R and excel. Turns out the problem was with excel and not with R
at all. That's a relief.

SOLVED




--
View this message in context: 
http://r.789695.n4.nabble.com/problem-with-factor-levels-tp4652006p4652019.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] computing marginal values based on multiple columns?

2012-12-04 Thread Simon
Hello all,

I have what feels like a simple problem, but I can't find an simple
answer. Consider this data frame:

 x - data.frame(sample1=c(35,176,182,193,124),
sample2=c(198,176,190,23,15), sample3=c(12,154,21,191,156),
class=c('a','a','c','b','c'))

 x
  sample1 sample2 sample3 class
1  35 198  12 a
2 176 176 154 a
3 182 190  21 c
4 193  23 191 b
5 124  15 156 c

Now I wish to know: for each sample, for values  20% of the sample mean,
what percentage of those are class a?

I want to end up with a table like:

   sample1 sample2 sample3
1  1.0 0 0.5

I can calculate this for an individual sample using this rather clumsy
expression:

length(which(x$sample1  mean(x$sample1)  x$class=='a')) /
length(which(x$sample1  mean(x$sample1)))

I'd normally propagate it across the data frame using apply, but I
can't because it depends on more than one column.

Any help much appreciated!

Cheers,

Simon

[[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] partial analisys of a time series

2012-12-04 Thread PIKAL Petr
Hi

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Antonio Silva
 Sent: Tuesday, December 04, 2012 10:26 AM
 To: R-help@r-project.org
 Subject: [R] partial analisys of a time series
 
 Dear list members
 
 I want to analyze separately the months of a time series. In other
 words, I want to plot and fit models for each month separately.
 
 Taking the example of
 http://a-little-book-of-r-for-time-
 series.readthedocs.org/en/latest/src/timeseries.html
 
 births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
 birthstimeseries - ts(births, frequency=12, start=c(1946,1))
 birthstimeseries
 plot.ts(birthstimeseries)
 birthstimeseriesHW - HoltWinters(birthstimeseries)
 plot(birthstimeseriesHW)
 
 How to proceed the plotting and HoltWinters smoothing considereing only
 Januarys, Februarys, etc. separately.

Split your data by months to a list, use lapply.

using zoo package

blist -split(birthstimeseries, months(as.Date(birthstimeseries)))
l.blist - lapply(blist, HoltWinters)

Regards
Petr





 
 Thanks in advance.
 
 Antonio Olinto
 
   [[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] problem with factor levels

2012-12-04 Thread PIKAL Petr
Hi

That is quite usual. Excel is so widespread that almost everybody assumes it 
shall not contain mistakes and behaves correctly. The contrary is true. 
Spreadsheet often guess what user have on mind and corrects values to fit 
such assumption, let aside mistakes in coded functions.

R expects it is used by clever and able people and performs just what you tell 
it to do, not more not less.

So whenever result does not fit your expectations, first proof your 
expectations.

Regards
Petr

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Jeremy.Shearman
 Sent: Tuesday, December 04, 2012 10:38 AM
 To: r-help@r-project.org
 Subject: Re: [R] problem with factor levels
 
 Oh, your skepticism was spot on!
 I was using excel to check the output (silly, but I am still in the
 process of moving from excel to R) and there was a discrepancy in the
 number of output from R and excel. Turns out the problem was with excel
 and not with R at all. That's a relief.
 
 SOLVED
 
 
 
 
 --
 View this message in context: http://r.789695.n4.nabble.com/problem-
 with-factor-levels-tp4652006p4652019.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] question about sum of (column) elements in R

2012-12-04 Thread Berend Hasselman

On 04-12-2012, at 08:59, T Bal wrote:

 Hi,
 I have the following data:
 
 0   12
 1   10
 1   4
 1   6
 1   7
 1   13
 2   21
 2   23
 2   20
 3   18
 3   17
 3   16
 3   27
 3   33
 4   11
 4   8
 4   19
 4   16
 4   9
 
 
 In this data file I would like to sum the numbers of second column which
 belong to the same number in the first column.
 So the output would be:
 
 0 12
 1 40
 2 64
3 111
etc.

Using Pascal's test object you can also use aggregate

aggregate(test$V2,by=list(test$V1),FUN=sum)

Berend

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 sum of (column) elements in R

2012-12-04 Thread Jose Iparraguirre
Hi,

Imagine the data you have is in a data frame, c, with columns a and b.
Then you can do this:

 res=tapply(b,c[,-2],sum) 
 res[is.na(res)]-0
 res
  0   1   2   3   4 
 12  40  64 111  63

Hope it helps,
José

José Iparraguirre
Chief Economist
Age UK



-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of T Bal
Sent: 04 December 2012 07:59
To: r-help@r-project.org
Subject: [R] question about sum of (column) elements in R

Hi,
I have the following data:

0   12
1   10
1   4
1   6
1   7
1   13
2   21
2   23
2   20
3   18
3   17
3   16
3   27
3   33
4   11
4   8
4   19
4   16
4   9


In this data file I would like to sum the numbers of second column which
belong to the same number in the first column.
So the output would be:

0 12
1 40
2 64
3 111
etc.

Thank you.

Kind regards,
T. Bal

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

A Star for Christmas

Kick start the festive season by attending one of Age UK’s Carol Concerts, A 
Star for Christmas. Taking place at Manchester Cathedral on Saturday 1 December 
and London’s St Pancras Church (opposite Euston Station) on Thursday 6 
December, they will feature special musical performances, readings by your 
favourite celebrities and carols, followed by mince pies and wine.
Tickets are priced at £20 full price/ £10 concessions. For more information, 
please visit http://www.ageuk.org.uk/astarforchristmas  or contact the 
Fundraising Events Team on
020 303 31725.

Age UK  Improving later life
www.ageuk.org.uk

 

---
Age UK is a registered charity and company limited by guarantee, (registered 
charity number 1128267, registered company number 6825798). 
Registered office: Tavis House, 1-6 Tavistock Square, London WC1H 9NA.

For the purposes of promoting Age UK Insurance, Age UK is an Appointed 
Representative of Age UK Enterprises Limited, Age UK is an Introducer 
Appointed Representative of JLT Benefit Solutions Limited and Simplyhealth 
Access for the purposes of introducing potential annuity and health 
cash plans customers respectively.  Age UK Enterprises Limited, JLT Benefit 
Solutions Limited and Simplyhealth Access are all authorised and 
regulated by the Financial Services Authority. 
--

This email and any files transmitted with it are confidential and intended 
solely for the use of the individual or entity to whom they are 
addressed. If you receive a message in error, please advise the sender and 
delete immediately.

Except where this email is sent in the usual course of our business, any 
opinions expressed in this email are those of the author and do not 
necessarily reflect the opinions of Age UK or its subsidiaries and associated 
companies. Age UK monitors all e-mail transmissions passing 
through its network and may block or modify mails which are deemed to be 
unsuitable.

Age Concern England (charity number 261794) and Help the Aged (charity number 
272786) and their trading and other associated companies merged 
on 1st April 2009.  Together they have formed the Age UK Group, dedicated to 
improving the lives of people in later life.  The three national 
Age Concerns in Scotland, Northern Ireland and Wales have also merged with Help 
the Aged in these nations to form three registered charities: 
Age Scotland, Age NI, Age Cymru.










__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] partial analisys of a time series

2012-12-04 Thread Antonio Silva
Thanks Petr

I thought there might be an equivalent for birthstimeseries[,1] if it were
a dataframe, but split function sounds great.

I could not reproduce the second line of your suggestion l.blist -
lapply(blist, HoltWinters). I receive the message: Error in
decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) : time
series has no or less than 2 periods

What could be going wrong?

Best regards

Antonio


2012/12/4 PIKAL Petr petr.pi...@precheza.cz

 Hi

  -Original Message-
  From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
  project.org] On Behalf Of Antonio Silva
  Sent: Tuesday, December 04, 2012 10:26 AM
  To: R-help@r-project.org
  Subject: [R] partial analisys of a time series
 
  Dear list members
 
  I want to analyze separately the months of a time series. In other
  words, I want to plot and fit models for each month separately.
 
  Taking the example of
  http://a-little-book-of-r-for-time-
  series.readthedocs.org/en/latest/src/timeseries.html
 
  births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
  birthstimeseries - ts(births, frequency=12, start=c(1946,1))
  birthstimeseries
  plot.ts(birthstimeseries)
  birthstimeseriesHW - HoltWinters(birthstimeseries)
  plot(birthstimeseriesHW)
 
  How to proceed the plotting and HoltWinters smoothing considereing only
  Januarys, Februarys, etc. separately.

 Split your data by months to a list, use lapply.

 using zoo package

 blist -split(birthstimeseries, months(as.Date(birthstimeseries)))
 l.blist - lapply(blist, HoltWinters)

 Regards
 Petr





 
  Thanks in advance.
 
  Antonio Olinto
 
[[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.




-- 
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[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] computing marginal values based on multiple columns?

2012-12-04 Thread Gerrit Eichner

Hello, Simon,

see below!


On Tue, 4 Dec 2012, Simon wrote:


Hello all,

I have what feels like a simple problem, but I can't find an simple
answer. Consider this data frame:


x - data.frame(sample1=c(35,176,182,193,124),

sample2=c(198,176,190,23,15), sample3=c(12,154,21,191,156),
class=c('a','a','c','b','c'))


x

 sample1 sample2 sample3 class
1  35 198  12 a
2 176 176 154 a
3 182 190  21 c
4 193  23 191 b
5 124  15 156 c

Now I wish to know: for each sample, for values  20% of the sample mean,
what percentage of those are class a?

I want to end up with a table like:

  sample1 sample2 sample3
1  1.0 0 0.5



I can't reproduce this result from your description above, but if I 
understand the latter correctly, maybe the following does what you want:


x.wo.class - subset( x, select = -class)
  # extract only the sample-columns

x.small.and.a - x.wo.class  0.2 * colMeans( x.wo.class)  x$class == a

apply( x.small.and.a, 2, function( xx) mean( x$class[ xx] == a))


 Hth  --  Gerrit



I can calculate this for an individual sample using this rather clumsy
expression:

length(which(x$sample1  mean(x$sample1)  x$class=='a')) /
length(which(x$sample1  mean(x$sample1)))

I'd normally propagate it across the data frame using apply, but I
can't because it depends on more than one column.

Any help much appreciated!

Cheers,

Simon

[[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] Resampling Help Needed

2012-12-04 Thread Suzen, Mehmet
You can use, 'sample' function for sampling  and may consider using
partition clustering for selecting your regions, see Cluster task view:

http://cran.r-project.org/web/views/Cluster.html

On 4 December 2012 00:53, KoopaTrooper ncoop...@tulane.edu wrote:
 I am using package ks() to build 3D representations of bird territories and
 calculate territory volume from spatial data (simply x, y, and z
 coordinates). What I want to do is determine at what sample size (#
 locations collected) does the territory volume stop increasing. This should
 give me an idea of the number of points needed for future seasons.

 So I have a couple of birds each with 200 spatial locations (x,y,z). I want
 to run the following code (see below), but have R calculate territory size
 100 times with 10 random points (no replacement), 100 times with 20 random
 points, 100 times with 30 random points, etc. I can figure out how to do
 this manually (i.e. create 100 individual files with 10 random points, 20
 random points, etc.) but I figure there must be a way to make my life
 easier. Any help would be appreciated. Even pointing me in the correct
 direction would be a big help. Thanks!

 Nathan

 #read data files (.csv's with 200 rows of x,y,z coordinates)
 a-read.csv(A_PW_ASY_M_LII_2011.csv)

 #calls the plug-in bandwidth estimator
 Ha - Hpi(a)

 #sets min/max grid size for each dimension
 minX-min(a$X)-25
 minY-min(a$Y)-25
 minZ-0

 maxX-max(a$X)+25
 maxY-max(a$Y)+25
 maxZ-max(a$Z)+5

 #creates kernel utilization distribution
 fhata - kde(x=a, H=Ha, binned=FALSE, xmin=c(minX,minY,minZ),
 xmax=c(maxX,maxY,maxZ))

 #calculates territory volume at 95% isopleth
 Vol95-contourSizes(fhata, cont=95)



 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Resampling-Help-Needed-tp4651973.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] Réponse automatique

2012-12-04 Thread j . boutet



Bonjour,



Je serais en congés jusqu'au Jeudi 6 Décembre.

Pour des raisons d'urgence, vous pourrez me contacter par téléphone au 06 46 34 
81 03.



Cordialement,



--



Jérôme Boutet

Conservatoire d'espaces naturels de Picardie

1, place Ginkgo - village Oasis

80 044 AMIENS cedex 

tél : 03 22 89 84 24 



 From: r-help-requ...@r-project.org
 Subject: R-help Digest, Vol 118, Issue 4
 Date: Tue, 04 Dec 2012 12:00:10 +0100

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

2012-12-04 Thread Michael Dewey

At 20:39 03/12/2012, Min Dong wrote:

Hi, I am a novice in R. It will be greatly appreciated if someone
can advise me with the following questions.


There are at least three packages available from CRAN (meta, metafor, 
rmeta) which draw forest plots so it would help us if you had told us 
which one you are using.




1) How to highlight reference range in forest plot? For example, if 1.5-2
is the reference range, I would like to have all the area between 1.5-2 to
be highlighted (such as in grey color).


Are you really using a forest plot for a meta-analysis or for some 
other purpose?



2) If I have handunds of objects, how to set the output plot to multiple
columns? For example, I have 500 objects included, and I want to have
objects 1-100 to the left (column #1), 101-200 in column#2 (next to
column#1)...etc, ect, how to do it?


It would be unusual to have a meta-analysis with so many studies 
included so I assume you are doing something else.



Thank you very much!

Mindy

[[alternative HTML version deleted]]


Michael Dewey
i...@aghmed.fsnet.co.uk
http://www.aghmed.fsnet.co.uk/home.html

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] partial analisys of a time series

2012-12-04 Thread R. Michael Weylandt
Try period.apply() from the xts package.

MW

On Tue, Dec 4, 2012 at 9:26 AM, Antonio Silva aolinto@gmail.com wrote:
 Dear list members

 I want to analyze separately the months of a time series. In other words, I
 want to plot and fit models for each month separately.

 Taking the example of
 http://a-little-book-of-r-for-time-series.readthedocs.org/en/latest/src/timeseries.html

 births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
 birthstimeseries - ts(births, frequency=12, start=c(1946,1))
 birthstimeseries
 plot.ts(birthstimeseries)
 birthstimeseriesHW - HoltWinters(birthstimeseries)
 plot(birthstimeseriesHW)

 How to proceed the plotting and HoltWinters smoothing considereing only
 Januarys, Februarys, etc. separately.

 Thanks in advance.

 Antonio Olinto

 [[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] partial analisys of a time series

2012-12-04 Thread PIKAL Petr
Hi

I am not an expert in time series. The problem is that resulting list is not 
time series any more. So you need to convert it again to time series and you 
need to give it frequency greater than 1.

something like
l.blist - lapply(blist, function (x) HoltWinters(ts(x, frequency=2)))

But this depends on how your data are structured after splitting them to months.

Regards
Petr


From: Antonio Silva [mailto:aolinto@gmail.com]
Sent: Tuesday, December 04, 2012 11:50 AM
To: PIKAL Petr
Cc: R-help@r-project.org
Subject: Re: [R] partial analisys of a time series

Thanks Petr

I thought there might be an equivalent for birthstimeseries[,1] if it were a 
dataframe, but split function sounds great.

I could not reproduce the second line of your suggestion l.blist - 
lapply(blist, HoltWinters). I receive the message: Error in 
decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) : time 
series has no or less than 2 periods

What could be going wrong?

Best regards

Antonio

2012/12/4 PIKAL Petr petr.pi...@precheza.czmailto:petr.pi...@precheza.cz
Hi

 -Original Message-
 From: r-help-boun...@r-project.orgmailto:r-help-boun...@r-project.org 
 [mailto:r-help-bounces@r-mailto:r-help-bounces@r-
 project.orghttp://project.org] On Behalf Of Antonio Silva
 Sent: Tuesday, December 04, 2012 10:26 AM
 To: R-help@r-project.orgmailto:R-help@r-project.org
 Subject: [R] partial analisys of a time series

 Dear list members

 I want to analyze separately the months of a time series. In other
 words, I want to plot and fit models for each month separately.

 Taking the example of
 http://a-little-book-of-r-for-time-
 series.readthedocs.org/en/latest/src/timeseries.htmlhttp://series.readthedocs.org/en/latest/src/timeseries.html

 births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
 birthstimeseries - ts(births, frequency=12, start=c(1946,1))
 birthstimeseries
 plot.ts(birthstimeseries)
 birthstimeseriesHW - HoltWinters(birthstimeseries)
 plot(birthstimeseriesHW)

 How to proceed the plotting and HoltWinters smoothing considereing only
 Januarys, Februarys, etc. separately.
Split your data by months to a list, use lapply.

using zoo package

blist -split(birthstimeseries, months(as.Date(birthstimeseries)))
l.blist - lapply(blist, HoltWinters)

Regards
Petr






 Thanks in advance.

 Antonio Olinto

   [[alternative HTML version deleted]]

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



--
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil


[[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] partial analisys of a time series

2012-12-04 Thread arun
HI,
You can subset by:
birthstimeseriesJan-subset(birthstimeseries,cycle(birthstimeseries)==1)
A.K.




- Original Message -
From: Antonio Silva aolinto@gmail.com
To: R-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 4:26 AM
Subject: [R] partial analisys of a time series

Dear list members

I want to analyze separately the months of a time series. In other words, I
want to plot and fit models for each month separately.

Taking the example of
http://a-little-book-of-r-for-time-series.readthedocs.org/en/latest/src/timeseries.html

births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
birthstimeseries - ts(births, frequency=12, start=c(1946,1))
birthstimeseries
plot.ts(birthstimeseries)
birthstimeseriesHW - HoltWinters(birthstimeseries)
plot(birthstimeseriesHW)

How to proceed the plotting and HoltWinters smoothing considereing only
Januarys, Februarys, etc. separately.

Thanks in advance.

Antonio Olinto

    [[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] Using multicores in R

2012-12-04 Thread moriah
Thanks for the help, 

Perhaps I should elaborate a bit, I am working on bioinformatics  project in
which I am trying to run a forward selection algorithm for machine learning
classification of two biological conditions. 
At each iteration I want to find the gene that in addition to those I have
found already does the best classification.

It looks something like this:

for (j in 1:5030)
  { 
  tp - 0;
  for (i in 1:5030)
  {
if (!(i %in% idx))
{

  classifier-naiveBayes(trn[,c(i,idx)], trn[,20118]) 
  tbl -table(predict(classifier, trn[,-20118]), trn[,20118])
  success - (tbl[[1]] +tbl[[4]])/(tbl[[1]] +tbl[[4]]+tbl[[2]]+tbl[[3]])

  if (success  tp)
  {
tp - success
ind - i
gene - names(trn)[i]
  }
}

  }
  idx - c(idx,ind)
  res - rbind(res, data.frame(Iteration=j,Success=tp*100,Gene=gene))
}

I am no expert when it comes to programming so I am not sure how can I
optimize my relatively primitive code in the best way...

Thanks,
Moriah





--
View this message in context: 
http://r.789695.n4.nabble.com/Using-multicores-in-R-tp4651808p4652034.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] do.call

2012-12-04 Thread Dominic Roye
Hello,

I have a problem with the do.call-function. I would like to merge
the values of more than 30 columns, but not in all of the rows exist
values, so with this commando i get a lot of ; or NA.

How get i only the merge of cells with a number?


datos$NEW - do.call(paste, c(datos[,19:53], sep = ;))


 $ NEW  : chr
218.0NA;;;NA;;;NA;NA;NA;NA;NA


I hope you can help me. Thanks!

Best regards,

Dominic

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 sum of (column) elements in R

2012-12-04 Thread arun
HI,
You can use either ?tapply(), ?aggregate(), ?ddply() from library(plyr)
dat1-read.table(text=
0  12
1  10
1  4
1  6
1  7
1  13
2  21
2  23
2  20
3  18
3  17
3  16
3  27
3  33
4  11
4  8
4  19
4  16
4  9
,sep=,header=FALSE)
 with(dat1,aggregate(dat1[,2],by=list(V1=dat1[,1]),sum))
 # V1   x
#1  0  12
#2  1  40
#3  2  64
#4  3 111
#5  4  63
A.K.




- Original Message -
From: T Bal studentt...@gmail.com
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 2:59 AM
Subject: [R] question about sum of (column) elements in R

Hi,
I have the following data:

0   12
1   10
1   4
1   6
1   7
1   13
2   21
2   23
2   20
3   18
3   17
3   16
3   27
3   33
4   11
4   8
4   19
4   16
4   9


In this data file I would like to sum the numbers of second column which
belong to the same number in the first column.
So the output would be:

0     12
1     40
2     64
        3     111
        etc.

Thank you.

Kind regards,
T. Bal

    [[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] Chi-squared test when observed near expected

2012-12-04 Thread David L Carlson
When you typed x as a command, R runs the command print(x). That function
produces a summary of the results which may include round off numbers to a
few decimal places to make them more readable. When you typed x$statistic,
you got the unrounded version of the result 5.6e-31 which I think you will
agree is pretty close to 0.

--
David L Carlson
Associate Professor of Anthropology
Texas AM University
College Station, TX 77843-4352

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Troy S
 Sent: Monday, December 03, 2012 3:41 PM
 To: r-help@r-project.org
 Subject: [R] Chi-squared test when observed near expected
 
 Dear UseRs,
 
 I'm running a chi-squared test where the expected matrix is the same as
 the
 observed, after rounding.  R reports a X-squared of zero with a p value
 of
 one.  I can justify this because any other result will deviate at least
 as
 much from the expected because what we observe is the expected, after
 rounding.  But the formula for X-squared, sum (O-E)^2/E gives a
 positive
 value.  What is the reason for X-Squared being zero in this case?
 
 Troy
 
  trial-as.table(matrix(c(26,16,13,7),ncol=2))
  x-chisq.test(trial)
  x
 
 
 
 data:  trial
 X-squared = 0, df = 1, p-value = 1
 
  x$expected
  A B
 A 26.41935 12.580645
 B 15.58065  7.419355
 
  x$statistic
X-squared
 5.596653e-31
  (x$observed-x$expected)^2/x$expected
 A   B
 A 0.006656426 0.013978495
 B 0.011286983 0.023702665
  sum((x$observed-x$expected)^2/x$expected)
 [1] 0.05562457
 
 
   [[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] julia language unfair comparisons against R

2012-12-04 Thread Suzen, Mehmet
Hello List,

Probably many of you aware of the Julia language
(http://julialang.org/), It is a promising project.

However it seems like R is very slow in their benchmarks. Very
important point they omit, they did not
use R's own JIT !  I had a feeling that R is mistreaded there :)
Also another important
point is that they all use for-loops in R instead of vectorized code!

Any thought on this?  Are there any other improvement one can do with
compiler package other then
'cmpfun'?

Best,
-m

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 make R randomForest model size smaller?

2012-12-04 Thread Liaw, Andy
Try the following:

set.seed(100)
rf1 - randomForest(Species ~ ., data=iris)
set.seed(100)
rf2 - randomForest(iris[1:4], iris$Species)
object.size(rf1)
object.size(rf2)
str(rf1)
str(rf2)

You can try it on your own data.  That should give you some hints about why the 
formula interface should be avoided with large datasets.

Andy

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of John Foreman
Sent: Monday, December 03, 2012 3:43 PM
To: r-help@r-project.org
Subject: [R] How do I make R randomForest model size smaller?

I've been training randomForest models on 7 million rows of data (41
features). Here's an example call:

myModel - randomForest(RESPONSE~., data=mydata, ntree=50, maxnodes=30)

I thought surely with only 50 trees and 30 terminal nodes that the memory
footprint of myModel would be small. But it's 65 megs in a dump file. The
object seems to be holding all sorts of predicted, actual, and vote data
from the training process.

What if I just want the forest and that's it? I want a tiny dump file that
I can load later to make predictions off of quickly. I feel like the forest
by itself shouldn't be all that large...

Anyone know how to strip this sucker down to just something I can make
predictions off of going forward?

[[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.
Notice:  This e-mail message, together with any attachme...{{dropped:11}}

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

2012-12-04 Thread R. Michael Weylandt
There been, that done.

http://stackoverflow.com/questions/9968578/speeding-up-julias-poorly-written-r-examples/10712158#10712158

MW

On Tue, Dec 4, 2012 at 2:34 PM, Suzen, Mehmet msu...@gmail.com wrote:
 Hello List,

 Probably many of you aware of the Julia language
 (http://julialang.org/), It is a promising project.

 However it seems like R is very slow in their benchmarks. Very
 important point they omit, they did not
 use R's own JIT !  I had a feeling that R is mistreaded there :)
 Also another important
 point is that they all use for-loops in R instead of vectorized code!

 Any thought on this?  Are there any other improvement one can do with
 compiler package other then
 'cmpfun'?

 Best,
 -m

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

2012-12-04 Thread anto.r
DT = data.frame(x=rep(c(a,b,c),each=3), y=c(1,3,6), v=1:9, w=3:11,
z=LETTERS[1:9]) 

If I understand you right, and you want to select all rows where v3 and
W10

with(DT, DT[which(v3  w10),])
  x y v w z
4 b 1 4 6 D
5 b 3 5 7 E
6 b 6 6 8 F
7 c 1 7 9 G

you can use colSums, rowSums  on this, but after removing the categorical
columns
Anto





--
View this message in context: 
http://r.789695.n4.nabble.com/Subsetting-columns-in-data-table-tp4652048p4652049.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] Winbugs from R

2012-12-04 Thread Günal Bilek
Hi,
I am trying to covert a Winbugs code into R code. Here is the winbugs code 
model{# model’s likelihoodfor (i in 1:n){time[i] ~ dnorm( mu[i], tau ) # 
stochastic componenent# link and linear predictormu[i] - beta0 + beta1 * 
cases[i] + beta2 * distance[i]}# prior distributionstau ~ dgamma( 0.01, 0.01 
)beta0 ~ dnorm( 0.0, 1.0E-4)beta1 ~ dnorm( 0.0, 1.0E-4)beta2 ~ dnorm( 0.0, 
1.0E-4)# definition of sigmas2-1/taus -sqrt(s2)# calculation of the sample 
variancefor (i in 1:n){ c.time[i]-time[i]-mean(time[]) }sy2 - inprod( 
c.time[], c.time[] )/(n-1)# calculation of Bayesian version R squaredR2B - 1 - 
s2/sy2# Expected y for a typical delivery timetypical.y - beta0 + beta1 * 
mean(cases[]) + beta2 * mean(distance[])}INITSlist( tau=1, beta0=1, beta1=0, 
beta2=0 )DATA (LIST)list( n=25,time = c(16.68, 11.5, 12.03, 14.88, 13.75, 
18.11, 8, 17.83,79.24, 21.5, 40.33, 21, 13.5, 19.75, 24, 29, 15.35,19, 9.5, 
35.1, 17.9, 52.32, 18.75, 19.83, 10.75),distance = c(560, 220, 340, 80, 150, 
330, 110, 210, 1460,605, 688, 215, 255, 462, 448, 776, 200, 132,36, 770, 140, 
810, 450, 635, 150),cases = c( 7, 3, 3, 4, 6, 7, 2, 7, 30, 5, 16, 10, 4, 6, 
9,10, 6, 7, 3, 17, 10, 26, 9, 8, 4) )

I want to do this in R. So, I copied the model and pasted into a txt file named 
reg. Here is the R code I used
time - 
c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.1,17.9,52.32,18.75,19.83,10.75)cases
 - c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)distance - 
c(560,220,340,80,150,330,110,210,1260,605,688,215,255,462,448,776,200,132,36,770,140,810,450,635,150)
data - list(time,cases,distance)
inits - function(){list(tau=1,beta1=0,beta2=0,beta3=0)}
sim - bugs(data, inits, model.file = 
C:/Users/Gunal/Desktop/dummy/reg.txt,parameters = c(beta1, 
beta2,beta3),n.chains = 3, n.iter = 1000,bugs.directory = 
D:/PROGRAMLAR/WinBUGS14/,debug=TRUE)

Winbugs is producing this error page
display(log)check(C:/Users/Gunal/Desktop/dummy/reg.txt)model is syntactically 
correctdata(C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/data.txt)data 
loadedcompile(3)variable n is not 
definedinits(1,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/inits1.txt)command 
#Bugs:inits cannot be executed (is greyed 
out)inits(2,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/inits2.txt)command 
#Bugs:inits cannot be executed (is greyed 
out)inits(3,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/inits3.txt)command 
#Bugs:inits cannot be executed (is greyed out)gen.inits()command 
#Bugs:gen.inits cannot be executed (is greyed 
out)thin.updater(1)update(500)command #Bugs:update cannot be executed (is 
greyed out)set(beta1)command #Bugs:set cannot be executed (is greyed 
out)set(beta2)command #Bugs:set cannot be executed (is greyed 
out)set(beta3)command #Bugs:set cannot be executed (is greyed 
out)set(deviance)command #Bugs:set cannot be executed (is greyed 
out)dic.set()command #Bugs:dic.set cannot be executed (is greyed 
out)update(500)command #Bugs:update cannot be executed (is greyed 
out)coda(*,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/coda)command #Bugs:coda 
cannot be executed (is greyed out)stats(*)command #Bugs:stats cannot be 
executed (is greyed out)dic.stats()
DIChistory(*,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/history.odc)command 
#Bugs:history cannot be executed (is greyed 
out)save(C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/log.odc)save(C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/log.txt)

Any help would be greatly appreciated. 
Cheers
Gunal 
[[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 for a function

2012-12-04 Thread anoumou
Hello all,
I need a help.
I am modeling a disease and a create a R function like that:

Lambda-function (x,date1,r,h,a){
  ndate1 - as.Date(date1, %d/%m/%Y)
  t1 - as.numeric(ndate1)
  x[order(x$i),]
  t -x[,t]
  i -x[,i]
  CONTAGIEUX -x[,CONTAGIEUX]
  while ( t1  min(t) ){
  for (i in 1:length(i) ){
{for (j in 1:CONTAGIEUX[length(CONTAGIEUX)]){
  res1[j] -(a*h)
  res2 -sum( res1[j])
}
 }
  lambda[i] - r*res2
  }
  }
x-data.frame(x,lambda)
x
}

on such data :

DATEi   Symptomes   t   Incubation  CONTAGIEUX
1   2009-04-29   Canada 13  14363   13  13
2   2009-05-01   Israel 2   14365   2   2
3   2009-05-09  argentina   1   14373   1   1
5   2009-05-09  australia   1   14373   1   1
6   2009-05-10  australia   1   14374   2   2
7   2009-04-29  Austria 1   14363   1   1
8   2009-04-30  Austria 1   14364   2   2
9   2009-05-01  Austria 1   14365   2   3
10  2009-05-02  Austria 1   14366   2   4
11  2009-05-03  Austria 1   14367   2   5
17  2009-05-09  Austria 1   14373   2   7
18  2009-05-10  Austria 1   14374   2   7
19  2009-05-08  brasil  4   14372   4   4
20  2009-05-09  brazil  6   14373   6   6
21  2009-05-10  brazil  6   14374   12  12
22  2009-05-02  canada  51  14366   51  51
23  2009-05-03  canada  85  14367   136 136
24  2009-05-04  canada  101 14368   186 237
31  2009-04-27  Canada  6   14361   6   6
32  2009-04-28  Canada  6   14362   6   6
33  2009-04-30  Canada  19  14364   25  25
34  2009-05-01  Canada  34  14365   53  59
35  2009-05-01  China,HongKong, SAR 1   14365   1   1
36  2009-05-02  China,HongKong, SAR 1   14366   2   2
37  2009-05-03  China,HongKong, SAR 1   14367   2   3
38  2009-05-04  China,HongKong, SAR 1   14368   2   4
44  2009-05-10  China,HongKong, SAR 1   14374   2   7
45  2009-05-04  colombia1   14368   1   1
46  2009-05-05  colombia1   14369   2   2
47  2009-05-06  colombia

But i do not get the results,i try by all means but i d'ont understant the
problem.
Thanks for your help.



--
View this message in context: 
http://r.789695.n4.nabble.com/Help-for-a-function-tp4652054.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] r function definition

2012-12-04 Thread FJ M

Are you using windows? If you are you may want to try to run your R code from a 
batch file:
REM on Microsoft Windows (adjust the path to R.exe as needed) 
C:\Program Files\R\R-2.13.2\bin\x64\R.exe CMD BATCH
C:\Users\Frank\Documents\R\Projects\Heinrich\Heinrich.txt 
C:\Users\Frank\Documents\R\Projects\Heinrich\Heinrich.out
PAUSE

Where ...Heinrich.txt is the following R code in a text file:
pdf(Heinrich.pdf)
## load Pearson package
library(PearsonDS)
## calculate probability density function Fig. 1
x - seq(0.0, +20.0, length=100)
hx- dpearsonIV(x,m=2.25,nu=5,location=17,scale=2)
k-0.0206631
plot(x, hx , type=l, lwd=2, tck=1, xlab=x value,  ylab=Density, main=J. 
Heinrich Fig. 1)
## calculate probability density function Fig 2
hx- dpearsonIV(x,m=0.75,nu=1.0,location=15,scale=0.5)
k-0.218991
plot(x, hx , type=l, lwd=2, tck=1, xlab=x value,  ylab=Density, main=J. 
Heinrich Fig. 2)
quit()
and Heinrich.out will contain the output of the results of the commands. Look 
at Heinrich.out for results and errors. Then correct the text file and rerun.
GL
Frank
Chicago

 Date: Mon, 3 Dec 2012 12:30:31 -0800
 From: q...@hotmail.com
 To: r-help@r-project.org
 Subject: [R] r function definition

 I am a very new R user. I am trying to write functons and debug functions.
 One problem for me is that I need to alwasy copy the whole function body and
 resubmit to R console every time I changed even one line of the function.
 Because I have long algorithm function, copying and pasting is very tedious
 for me. I assume if I save the function files, R should be able to just use
 the new function body since it is a scripting language. Can somebody let me
 know his best practice of using R function?



 --
 View this message in context: 
 http://r.789695.n4.nabble.com/r-function-definition-tp4651943.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] help with if statement

2012-12-04 Thread anto.r
DT-data.frame(time=c(0,1,5,24,36,48,72),DV=seq(0,60,10))
  time DV
10  0
21 10
35 20
4   48 30
5   84 40
6   96 50
7  120 60
You want to add 24 to values that are =24 in 'time'

DT[DT$time=24,'time']-DT[DT$time=24,'time']+24
  time DV
10  0
21 10
35 20
4   48 30
5   60 40
6   72 50
7   96 60

Is this what you were looking for?
A.




--
View this message in context: 
http://r.789695.n4.nabble.com/help-with-if-statement-tp4652015p4652055.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] partial analisys of a time series

2012-12-04 Thread arun
HI,
I am getting an error message:
l.blist-lapply(blist,HoltWinters)
#Error in decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) 
: 
 # time series has no or less than 2 periods
A.K.



- Original Message -
From: PIKAL Petr petr.pi...@precheza.cz
To: Antonio Silva aolinto@gmail.com; R-help@r-project.org 
R-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 5:07 AM
Subject: Re: [R] partial analisys of a time series

Hi

 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
 project.org] On Behalf Of Antonio Silva
 Sent: Tuesday, December 04, 2012 10:26 AM
 To: R-help@r-project.org
 Subject: [R] partial analisys of a time series
 
 Dear list members
 
 I want to analyze separately the months of a time series. In other
 words, I want to plot and fit models for each month separately.
 
 Taking the example of
 http://a-little-book-of-r-for-time-
 series.readthedocs.org/en/latest/src/timeseries.html
 
 births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
 birthstimeseries - ts(births, frequency=12, start=c(1946,1))
 birthstimeseries
 plot.ts(birthstimeseries)
 birthstimeseriesHW - HoltWinters(birthstimeseries)
 plot(birthstimeseriesHW)
 
 How to proceed the plotting and HoltWinters smoothing considereing only
 Januarys, Februarys, etc. separately.

Split your data by months to a list, use lapply.

using zoo package

blist -split(birthstimeseries, months(as.Date(birthstimeseries)))
l.blist - lapply(blist, HoltWinters)

Regards
Petr





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

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


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


Re: [R] do.call

2012-12-04 Thread arun
Hi,
Try this:
set.seed(15)
datos-as.data.frame(matrix(sample(c(1:20,NA),30,replace=TRUE),ncol=6))
 do.call(paste,c(na.omit(datos),sep=;))
#[1] 5;18;14;10;17;3  14;15;15;3;2;20  8;18;19;17;12;11

A.K.



- Original Message -
From: Dominic Roye dominic.r...@gmail.com
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 7:38 AM
Subject: [R] do.call

Hello,

I have a problem with the do.call-function. I would like to merge
the values of more than 30 columns, but not in all of the rows exist
values, so with this commando i get a lot of ; or NA.

How get i only the merge of cells with a number?


datos$NEW - do.call(paste, c(datos[,19:53], sep = ;))


$ NEW                              : chr
218.0NA;;;NA;;;NA;NA;NA;NA;NA


I hope you can help me. Thanks!

Best regards,

Dominic

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] partial analisys of a time series

2012-12-04 Thread arun
Hi,
If the frequency is 1, the error message will  be gone.
For e.g.
birthstimeseriesJanFeb-subset(birthstimeseries,cycle(birthstimeseries)==c(1,2))
 birthstimeseriesJanFeb1-ts(birthstimeseriesJanFeb,frequency=2,start=c(1946,1))
 plot.ts(birthstimeseriesJanFeb1)
 birthstimeseriesJanFebHW-HoltWinters(birthstimeseriesJanFeb1)
 plot(birthstimeseriesJanFebHW)
A.K.




- Original Message -
From: Antonio Silva aolinto@gmail.com
To: PIKAL Petr petr.pi...@precheza.cz
Cc: R-help@r-project.org R-help@r-project.org
Sent: Tuesday, December 4, 2012 5:50 AM
Subject: Re: [R] partial analisys of a time series

Thanks Petr

I thought there might be an equivalent for birthstimeseries[,1] if it were
a dataframe, but split function sounds great.

I could not reproduce the second line of your suggestion l.blist -
lapply(blist, HoltWinters). I receive the message: Error in
decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) : time
series has no or less than 2 periods

What could be going wrong?

Best regards

Antonio


2012/12/4 PIKAL Petr petr.pi...@precheza.cz

 Hi

  -Original Message-
  From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
  project.org] On Behalf Of Antonio Silva
  Sent: Tuesday, December 04, 2012 10:26 AM
  To: R-help@r-project.org
  Subject: [R] partial analisys of a time series
 
  Dear list members
 
  I want to analyze separately the months of a time series. In other
  words, I want to plot and fit models for each month separately.
 
  Taking the example of
  http://a-little-book-of-r-for-time-
  series.readthedocs.org/en/latest/src/timeseries.html
 
  births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
  birthstimeseries - ts(births, frequency=12, start=c(1946,1))
  birthstimeseries
  plot.ts(birthstimeseries)
  birthstimeseriesHW - HoltWinters(birthstimeseries)
  plot(birthstimeseriesHW)
 
  How to proceed the plotting and HoltWinters smoothing considereing only
  Januarys, Februarys, etc. separately.

 Split your data by months to a list, use lapply.

 using zoo package

 blist -split(birthstimeseries, months(as.Date(birthstimeseries)))
 l.blist - lapply(blist, HoltWinters)

 Regards
 Petr





 
  Thanks in advance.
 
  Antonio Olinto
 
        [[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.




-- 
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

    [[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] Different results from random.Forest with test option and using predict function

2012-12-04 Thread Liaw, Andy
Without data to reproduce what you saw, we can only guess.

One possibility is due to tie-breaking.  There are several places where ties 
can occur and are broken at random, including at the prediction step.  One 
difference between the two ways of doing prediction is that when it's all done 
within randomForest(), the test set prediction is performed as each tree is 
grown.  If there is any tie that needs to be broken at any prediction step, it 
will affect the RNG stream used by the subsequent tree growing step.

You can also inspect/compare the forest components of the randomForest 
objects to see if they are the same.  At least the first tree in both should be 
identical.

Andy

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
Behalf Of tdbuskirk
Sent: Monday, December 03, 2012 6:31 PM
To: r-help@r-project.org
Subject: [R] Different results from random.Forest with test option and using 
predict function

Hello R Gurus,

I am perplexed by the different results I obtained when I ran code like
this:
set.seed(100)
test1-randomForest(BinaryY~., data=Xvars, trees=51, mtry=5, seed=200)
predict(test1, newdata=cbind(NewBinaryY, NewXs), type=response)

and this code:
set.seed(100)
test2-randomForest(BinaryY~., data=Xvars, trees=51, mtry=5, seed=200,
xtest=NewXs, ytest=NewBinarY)

The confusion matrices for the two forests I thought would be the same by
virtue of the same seed settings, but they differ as do the predicted values
as well as the votes.  At first I thought it was just the way ties were
broken, so I changed the number of trees to an odd number so there are no
ties anymore.  

Can anyone shed light on what I am hoping is a simple oversight?  I just
can't figure out why the results of the predictions from these two forests
applied to the NewBinaryYs and NewX data sets would not be the same.

Thanks for any hints and help.

Sincerely,

Trent Buskirk



--
View this message in context: 
http://r.789695.n4.nabble.com/Different-results-from-random-Forest-with-test-option-and-using-predict-function-tp4651970.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.
Notice:  This e-mail message, together with any attachme...{{dropped:11}}

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

2012-12-04 Thread YAddo
Thanks, Rui and David!



--
View this message in context: 
http://r.789695.n4.nabble.com/Histogram-plot-help-tp4651958p4652065.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] Spatial AND Temporal Dependency with lme, or other function?

2012-12-04 Thread Anne-Sophie Julien
Hi,

I've got a new dataset which I don't know how to analyze with R. My knowledge 
about R is limited for this kind of problem. I've tried to find a solution with 
some spatio temporel packages and lme/lmer functions, but didn't find any 
similar example.

I've got 10 locations on a coast on which the number of walrus was noted. I've 
got a distance measure between each location. We also have an indicator of the 
presence of tourists, and the number of seals. These information on each 
location were noted 10 times ( 2 summers * 5 weeks). We want to know how the 
number of walrus is affected by the tourists, the number of seals, the week of 
the summer and the location. For example, we think that they are moving east at 
the end of the summer, when there are more tourist in the west, and less seals 
in the east.

My problem is how to take into account the spatial and temporal dependency. The 
locations are beside one another, so the number of walrus in one location is 
closely related to the closest location. I would like to take into account the 
distance between each site for the dependency. And we took 10 measures in time 
at each location, so I should take into account this correlation too.

I've looked at lme function to model the dependency. It's correct with only 
spatial or temporal, but I don't know how to take into account both 
dependencies at the same time.

If you know how to do it with lme, or with any other function/package, let me 
know!

Thanks in advance, and have a nice day!

Sophie


[[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] model selection with spg and AIC (or, convert list to fitted model object)

2012-12-04 Thread Ravi Varadhan
Adam,

Getting the variance of MLE estimator when the true parameter is on the 
boundary is a very difficult problem.  It is known that the standard bootstrap 
does not work.  There are some sub-sampling approaches (Springer book:  
Politis, Romano, Wolff), but I am not an expert on this.

However, an important question is how do we know that the truth is supposedly 
on the boundary.  If you can assume that the truth is in the interior, then you 
can use standard approaches: observed hessian or bootstrap to get standard 
errors.

Best,
Ravi

From: Adam Zeilinger [mailto:ad...@ucr.edu]
Sent: Sunday, December 02, 2012 5:04 PM
To: Ravi Varadhan
Cc: Adam Zeilinger (zeil0...@umn.edu); r-help@r-project.org
Subject: Re: [R] model selection with spg and AIC (or, convert list to fitted 
model object)

Dear Ravi,

Thank you so much for the help.  I switched to using the optimx function but I 
continue to use the spg method (for the most part) because I found that only 
spg consistently converges give different datasets.  I also decided to use AIC 
rather that a likelihood ratio test.

I have a new question.  I would like to construct 95% confidence intervals for 
the parameter estimates from the best model.  From a previous R Help thread, 
you said that it was a bad idea to use the Hessian matrix computed from 
optim/optimx or hessian() [numDeriv package] when the optimization is 
constrained and parameter estimates are on the boundary because the MLE is 
likely at a local minimum:

http://tolstoy.newcastle.edu.au/R/e15/help/11/09/6673.html

In the same thread, you suggest using the Hessian matrix from augmented 
Lagrangian optimization with auglag() [alabama package] (with some caveats).  I 
would like to construct 95% CI with auglag, but I don't understand how to write 
the inequality constraints (hin) function.  Could you please help me write the 
hin function?

Below is R code for my MLE problem, using data that results in parameter 
estimates on the boundary, and my unsuccessful attempt at auglag() 
optimization.  Note: I have gradient functions for NLL1 and NLL2 but they're 
very large and don't seem to improve optimization, so they are not included 
here.  I can supply them if needed for the auglag() function.  Also, I am 
running R 2.15.2 on a Windows 7 64-bit machine.

##

library(optimx)
library(alabama)

# define multinomial distribution
dmnom2 - function(x,prob,log=FALSE) {
  r - lgamma(sum(x) + 1) + sum(x * log(prob) - lgamma(x + 1))
  if (log) r else exp(r)
}

# data frame y


y - structure(list(t = c(0.167, 0.5, 1, 12, 18, 24, 36), n1 = c(1L,

1L, 1L, 8L, 9L, 10L, 12L), n2 = c(0L, 1L, 2L, 6L, 5L, 3L, 2L),

n3 = c(13L, 12L, 11L, 0L, 0L, 1L, 0L)), .Names = c(t, n1,

n2, n3), class = data.frame, row.names = 36:42)

# Negative log-likelihood functions
NLL1 - function(par, y) {
  p1 - par[1]
  p2 - p1
  mu1 - par[2]
  mu2 - mu1
  t - y$t
  n1 - y$n1
  n2 - y$n2
  n3 - y$n3
  P1 - (p1*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*((-mu2)*(mu2 - p1 + p2) +
mu1*(mu2 + 2*p2)) - mu2*sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))) -
exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
mu2*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2*t)*mu2*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)/
exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2
  P2 - (p2*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*(-mu1^2 + 2*mu2*p1 +
mu1*(mu2 - p1 + p2)) - mu1*sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2))) -
exp(sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)))*t)*
mu1*sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2))) +
2*exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2*t)*mu1*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2)/
exp((1/2)*(mu1 + mu2 + p1 + p2 + sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2*t)/(2*(mu2*p1 + mu1*(mu2 + p2))*
sqrt((mu1 + mu2 + p1 + p2)^2 - 4*(mu2*p1 + mu1*(mu2 + p2
  P3 - 1 - P1 - P2
  p.all - c(P1, P2, P3)
  if(all(p.all  0  p.all  1)) -sum(dmnom2(c(n1, n2, n3), prob = p.all, log = 
TRUE)) else 1e07
}

NLL2 - function(par, y) {
  p1 - par[1]
  p2 - par[2]
  mu1 - par[3]
  mu2 - par[4]
  t - y$t
  n1 - y$n1
  n2 - y$n2
  n3 - y$n3
  P1 - (p1*((-1 + exp(sqrt((mu1 + mu2 + p1 + p2)^2 -
4*(mu2*p1 + mu1*(mu2 + p2)))*t))*((-mu2)*(mu2 - p1 + p2) +
mu1*(mu2 + 2*p2)) - mu2*sqrt((mu1 + mu2 + p1 + 

[R] lme: subject-specific slopes.

2012-12-04 Thread John Sorkin
I am running a random intercept random slope regression:
 
fitRIRT - lme(echogen~time,random=~  
1+time|subject,data=repeatdata,na.action=na.omit)
summary(fitRIRT)
 
I would like to get the subject-specific slopes, i.e. the slope that the model 
computes for each subject. If I have 10-subjects I should have 10-slopes. I 
don't see the slope when I look at the items listed in 
names(summary(fitRIRT)
nor when I look at the items listed in 
names(fitRIRT)
 
Thanks 
John
 
 

 
 
 
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)
Confidentiality Statement:
This email message, including any attachments, is for the sole use of the 
intended recipient(s) and may contain confidential and privileged information.  
Any unauthorized use, disclosure or distribution is prohibited.  If you are not 
the intended recipient, please contact the sender by reply email and destroy 
all copies of the original message. 
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lme: subject-specific slopes.

2012-12-04 Thread John Sorkin
Ken,
Thank you for your help. ranef(fitRIRT) does not give me what I expect. The 
subject-specific slopes, and subject-specific intercepts are not anywhere close 
to what I would expect them to be; the mean of the subject-specfic values 
should be close to those reported by 
summary(fitRIRT) and they are not. As you will see by examining the material 
below, the subject-specific slopes are off by many order of magnitude. The 
intercepts are also far from the value reported in summary(fitRIRT). Do you 
have any thoughts?
Thanks,
John
 
 fitRIRT - lme(echogen~time,random=~  
 1+time|subject,data=repeatdata,na.action=na.omit) summary(fitRIRT)Linear 
 mixed-effects model fit by REML
 Data: repeatdata 
  AIC  BIClogLik
  495.097 507.2491 -241.5485

Random effects:
 Formula: ~1 + time | subject
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr  
(Intercept) 1.917511e+01 (Intr)
time2.032276e-04 0 
Residual1.044601e+01   

Fixed effects: echogen ~ time 
   Value Std.Error DF   t-value p-value
(Intercept) 64.54864  4.258235 32 15.158543  0.
time 0.35795  0.227080 32  1.576307  0.1248
 Correlation: 
 (Intr)
time -0.242

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max 
-1.61362755 -0.52710871  0.02948008  0.41793322  1.77340082 

Number of Observations: 58
Number of Groups: 25  ranef(fitRIRT)   (Intercept)  time
1-3.278112  2.221016e-09
2   -35.400618  4.314995e-08
311.493110 -6.797543e-09
4   -16.209586 -7.070834e-08
5 3.585227 -2.389705e-08
6 1.614320 -1.967700e-09
7 8.346905  5.827094e-08
830.917812 -3.768584e-08
9-0.394101 -9.158251e-09
104.437509 -4.057971e-08
11   31.956597 -2.126275e-08
12   41.567402 -4.853942e-08
13  -10.723993  1.307152e-08
14   -4.554837  5.551908e-09
15   -4.554501  4.815086e-08
16   13.296985 -3.743967e-08
17   -8.255439  1.733238e-08
18  -21.317239  2.203885e-08
19  -13.480159  2.194016e-08
20  -13.044766  2.269168e-08
21   11.639198 -1.418706e-08
22  -27.457388 -1.154099e-08
232.194001 -5.509119e-09
24   -3.992646  7.682188e-08
251.614320 -1.967700e-09
 
 


 
 
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing) Kenneth 
Frost kfr...@wisc.edu 12/4/2012 10:44 AM 
I think this might be close to what you want. ranef(fitRIRT)

On 12/04/12, John Sorkin  wrote:
 I am running a random intercept random slope regression:
 
 fitRIRT - lme(echogen~time,random=~ 
 1+time|subject,data=repeatdata,na.action=na.omit)
 summary(fitRIRT)
 
 I would like to get the subject-specific slopes, i.e. the slope that the 
 model computes for each subject. If I have 10-subjects I should have 
 10-slopes. I don't see the slope when I look at the items listed in 
 names(summary(fitRIRT)
 nor when I look at the items listed in 
 names(fitRIRT)
 
 Thanks 
 John
 
 
 
 
 
 
 John David Sorkin M.D., Ph.D.
 Chief, Biostatistics and Informatics
 University of Maryland School of Medicine Division of Gerontology
 Baltimore VA Medical Center
 10 North Greene Street
 GRECC (BT/18/GR)
 Baltimore, MD 21201-1524
 (Phone) 410-605-7119
 (Fax) 410-605-7913 (Please call phone number above prior to faxing)
 Confidentiality Statement:
 This email message, including any attachments, is for the sole use of the 
 intended recipient(s) and may contain confidential and privileged 
 information. Any unauthorized use, disclosure or distribution is prohibited. 
 If you are not the intended recipient, please contact the sender by reply 
 email and destroy all copies of the original message. 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

--
Kenneth Frost
Graduate Research Assistant - Dept. of Plant Pathology
University of Wisconsin - Madison
Lab: (608) 262-9914
Mobile: (608) 556-9637
kfr...@wisc.edu


Confidentiality Statement:
This email message, including any attachments, is for the sole use of the 
intended recipient(s) and may contain confidential and privileged information.  
Any unauthorized use, disclosure or distribution is prohibited.  If you are not 
the intended recipient, please contact the sender by reply email and destroy 
all copies of the original message. 
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] lme: subject-specific slopes.

2012-12-04 Thread Kenneth Frost
I think the random effects represent the subject adjustments to the population 
averages. You may have to do the addition yourself to get the subject specific 
slopes and intercepts. Someone will hopefully correct me if I'm wrong.

On 12/04/12, John Sorkin  wrote:
 
 
 
 Ken,
 Thank you for your help. ranef(fitRIRT) does not give me what I expect. The 
 subject-specific slopes, and subject-specific intercepts are not anywhere 
 close to what I would expect them to be; the mean of the subject-specfic 
 values should be close to those reported by 
 summary(fitRIRT) and they are not. As you will see by examining the material 
 below, the subject-specific slopes are off by many order of magnitude. The 
 intercepts are also far from the value reported in summary(fitRIRT). Do you 
 have any thoughts?
 Thanks,
 John
 
 
  fitRIRT - lme(echogen~time,random=~ 
  1+time|subject,data=repeatdata,na.action=na.omit)  summary(fitRIRT) Linear 
  mixed-effects model fit by REML Data: repeatdata AIC BIC logLik 495.097 
  507.2491 -241.5485 Random effects: Formula: ~1 + time | subject Structure: 
  General positive-definite, Log-Cholesky parametrization StdDev Corr 
  (Intercept) 1.917511e+01 (Intr) time 2.032276e-04 0 Residual 1.044601e+01 
  Fixed effects: echogen ~ time Value Std.Error DF t-value p-value 
  (Intercept) 64.54864 4.258235 32 15.158543 0. time 0.35795 0.227080 32 
  1.576307 0.1248 Correlation: (Intr) time -0.242 Standardized Within-Group 
  Residuals: Min Q1 Med Q3 Max -1.61362755 -0.52710871 0.02948008 0.41793322 
  1.77340082 Number of Observations: 58 Number of Groups: 25  ranef(fitRIRT) 
  (Intercept) time 1 -3.278112 2.221016e-09 2 -35.400618 4.314995e-08 3 
  11.493110 -6.797543e-09 4 -16.209586 -7.070834e-08 5 3.585227 -2.389705e-08 
  6 1.614320 -1.967700e-09 7 8.346905 5.827094e-08 8 30.917812 -3.768584e-08 
  9 -0.394101 -9.158251e-09 10 4.437509 -4.057971e-08 11 31.956597 
  -2.126275e-08 12 41.567402 -4.853942e-08 13 -10.723993 1.307152e-08 14 
  -4.554837 5.551908e-09 15 -4.554501 4.815086e-08 16 13.296985 -3.743967e-08 
  17 -8.255439 1.733238e-08 18 -21.317239 2.203885e-08 19 -13.480159 
  2.194016e-08 20 -13.044766 2.269168e-08 21 11.639198 -1.418706e-08 22 
  -27.457388 -1.154099e-08 23 2.194001 -5.509119e-09 24 -3.992646 
  7.682188e-08 25 1.614320 -1.967700e-09
 
 
 
 
 
 
 
 John David Sorkin M.D., Ph.D.
 Chief, Biostatistics and Informatics
 University of Maryland School of Medicine Division of Gerontology
 Baltimore VA Medical Center
 10 North Greene Street
 GRECC (BT/18/GR)
 Baltimore, MD 21201-1524
 (Phone) 410-605-7119
 (Fax) 410-605-7913 (Please call phone number above prior to faxing) 
 Kenneth Frost kfr...@wisc.edu 12/4/2012 10:44 AM 
 I think this might be close to what you want. ranef(fitRIRT)
 
 On 12/04/12, John Sorkin wrote:
  I am running a random intercept random slope regression:
  
  fitRIRT - lme(echogen~time,random=~ 
  1+time|subject,data=repeatdata,na.action=na.omit)
  summary(fitRIRT)
  
  I would like to get the subject-specific slopes, i.e. the slope that the 
  model computes for each subject. If I have 10-subjects I should have 
  10-slopes. I don't see the slope when I look at the items listed in 
  names(summary(fitRIRT)
  nor when I look at the items listed in 
  names(fitRIRT)
  
  Thanks 
  John
  
  
  
  
  
  
  John David Sorkin M.D., Ph.D.
  Chief, Biostatistics and Informatics
  University of Maryland School of Medicine Division of Gerontology
  Baltimore VA Medical Center
  10 North Greene Street
  GRECC (BT/18/GR)
  Baltimore, MD 21201-1524
  (Phone) 410-605-7119
  (Fax) 410-605-7913 (Please call phone number above prior to faxing)
  Confidentiality Statement:
  This email message, including any attachments, is for the sole use of the 
  intended recipient(s) and may contain confidential and privileged 
  information. Any unauthorized use, disclosure or distribution is 
  prohibited. If you are not the intended recipient, please contact the 
  sender by reply email and destroy all copies of the original message. 
  
  
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 --
 Kenneth Frost
 Graduate Research Assistant - Dept. of Plant Pathology
 University of Wisconsin - Madison
 Lab: (608) 262-9914
 Mobile: (608) 556-9637
 kfr...@wisc.edu
 
 
 
 
 Confidentiality Statement: 
 
 This email message, including any attachments, is for the sole use of the 
 intended recipient(s) and may contain confidential and privileged 
 information. Any unauthorized use, disclosure or distribution is prohibited. 
 If you are not the intended recipient, please contact the sender by reply 
 email and destroy all copies of the original message.

--
Kenneth Frost
Graduate Research Assistant - Dept. of Plant Pathology
University of Wisconsin 

Re: [R] lme: subject-specific slopes.

2012-12-04 Thread John Sorkin
Yes, you are correct. 
Thanks,
John

 
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing) Kenneth 
Frost kfr...@wisc.edu 12/4/2012 11:07 AM 
I think the random effects represent the subject adjustments to the population 
averages. You may have to do the addition yourself to get the subject specific 
slopes and intercepts. Someone will hopefully correct me if I'm wrong.

On 12/04/12, John Sorkin  wrote:
 
 
 
 Ken,
 Thank you for your help. ranef(fitRIRT) does not give me what I expect. The 
 subject-specific slopes, and subject-specific intercepts are not anywhere 
 close to what I would expect them to be; the mean of the subject-specfic 
 values should be close to those reported by 
 summary(fitRIRT) and they are not. As you will see by examining the material 
 below, the subject-specific slopes are off by many order of magnitude. The 
 intercepts are also far from the value reported in summary(fitRIRT). Do you 
 have any thoughts?
 Thanks,
 John
 
 
  fitRIRT - lme(echogen~time,random=~ 
  1+time|subject,data=repeatdata,na.action=na.omit)  summary(fitRIRT) Linear 
  mixed-effects model fit by REML Data: repeatdata AIC BIC logLik 495.097 
  507.2491 -241.5485 Random effects: Formula: ~1 + time | subject Structure: 
  General positive-definite, Log-Cholesky parametrization StdDev Corr 
  (Intercept) 1.917511e+01 (Intr) time 2.032276e-04 0 Residual 1.044601e+01 
  Fixed effects: echogen ~ time Value Std.Error DF t-value p-value 
  (Intercept) 64.54864 4.258235 32 15.158543 0. time 0.35795 0.227080 32 
  1.576307 0.1248 Correlation: (Intr) time -0.242 Standardized Within-Group 
  Residuals: Min Q1 Med Q3 Max -1.61362755 -0.52710871 0.02948008 0.41793322 
  1.77340082 Number of Observations: 58 Number of Groups: 25  ranef(fitRIRT) 
  (Intercept) time 1 -3.278112 2.221016e-09 2 -35.400618 4.314995e-08 3 
  11.493110 -6.797543e-09 4 -16.209586 -7.070834e-08 5 3.585227 -2.389705e-08 
  6 1.614320 -1.967700e-09 7 8.346905 5.827094e-08 8 30.917812 -3.768584e-08 
  9 -0.394101 -9.158251e-09 10 4.437509 -4.057971e-08 11 31.956597 
  -2.126275e-08 12 41.567402 -4.853942e-08 13 -10.723993 1.307152e-08 14 
  -4.554837 5.551908e-09 15 -4.554501 4.815086e-08 16 13.296985 -3.743967e-08 
  17 -8.255439 1.733238e-08 18 -21.317239 2.203885e-08 19 -13.480159 
  2.194016e-08 20 -13.044766 2.269168e-08 21 11.639198 -1.418706e-08 22 
  -27.457388 -1.154099e-08 23 2.194001 -5.509119e-09 24 -3.992646 
  7.682188e-08 25 1.614320 -1.967700e-09
 
 
 
 
 
 
 
 John David Sorkin M.D., Ph.D.
 Chief, Biostatistics and Informatics
 University of Maryland School of Medicine Division of Gerontology
 Baltimore VA Medical Center
 10 North Greene Street
 GRECC (BT/18/GR)
 Baltimore, MD 21201-1524
 (Phone) 410-605-7119
 (Fax) 410-605-7913 (Please call phone number above prior to faxing) 
 Kenneth Frost kfr...@wisc.edu 12/4/2012 10:44 AM 
 I think this might be close to what you want. ranef(fitRIRT)
 
 On 12/04/12, John Sorkin wrote:
  I am running a random intercept random slope regression:
  
  fitRIRT - lme(echogen~time,random=~ 
  1+time|subject,data=repeatdata,na.action=na.omit)
  summary(fitRIRT)
  
  I would like to get the subject-specific slopes, i.e. the slope that the 
  model computes for each subject. If I have 10-subjects I should have 
  10-slopes. I don't see the slope when I look at the items listed in 
  names(summary(fitRIRT)
  nor when I look at the items listed in 
  names(fitRIRT)
  
  Thanks 
  John
  
  
  
  
  
  
  John David Sorkin M.D., Ph.D.
  Chief, Biostatistics and Informatics
  University of Maryland School of Medicine Division of Gerontology
  Baltimore VA Medical Center
  10 North Greene Street
  GRECC (BT/18/GR)
  Baltimore, MD 21201-1524
  (Phone) 410-605-7119
  (Fax) 410-605-7913 (Please call phone number above prior to faxing)
  Confidentiality Statement:
  This email message, including any attachments, is for the sole use of the 
  intended recipient(s) and may contain confidential and privileged 
  information. Any unauthorized use, disclosure or distribution is 
  prohibited. If you are not the intended recipient, please contact the 
  sender by reply email and destroy all copies of the original message. 
  
  
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 --
 Kenneth Frost
 Graduate Research Assistant - Dept. of Plant Pathology
 University of Wisconsin - Madison
 Lab: (608) 262-9914
 Mobile: (608) 556-9637
 kfr...@wisc.edu
 
 
 
 
 Confidentiality Statement: 
 
 This email message, including any 

[R] How to reset R settings to original state other than remove .Rdata and .Rhistory

2012-12-04 Thread Tian3507

Dear all,
Do you know how to return all the R settings to original state? Other than 
.Rdata and .Rhistory
 
Some weid thing happened to my machine.  I was trying to get shaded confidence 
band ploted using survplot from rms liberary.
 
It worked on one machine, but not on the other. I tried unstall R and reinstall 
R and remove. Rdata and .Rhisotory. Nothing helped so far.
 
Thanks for your input.
Best regards, 
Hong
 
Codes are
library(rms)
data generation
n - 1000
set.seed(731)
age - 50 + 12*rnorm(n)
label(age) - Age
sex - factor(sample(c('male','female'), n, TRUE))
cens - 15*runif(n)
h - .02*exp(.04*(age-50)+.8*(sex=='female'))
dt - -log(runif(n))/h
label(dt) - 'Follow-up Time'
e - ifelse(dt = cens,1,0)
dt - pmin(dt, cens)
units(dt) - 'Year'
dd - datadist(age, sex)
options(datadist='dd')
S - Surv(dt,e)
tte-data.frame(dt,e,age,sex)
male group KM curve with bands#
f - survfit(Surv(dt,e)~1, data=tte,subset=sex=='male')
survplot(f,conf=bands)
female group KM curve with bands#
g - survfit(Surv(dt,e)~1, data=tte,subset=sex!='male')
survplot(g,conf=bands, col=3, add=T)  
  
[[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] Solve system of equations (nleqslv) only returns origin

2012-12-04 Thread Alicia Ellis
I'm solving 4 complex equations simultaneously.  Code is below.  The code
returns only zero's for the solution though there should also be a non-zero
result.  I'm pretty confident that the equations are correct because they
are straight from a published paper and I checked them pretty thoroughly.
The parameter values I used are from the published paper as well.  Any
suggestions for how I can find the maximum non-zero solution instead of
going straight to the minimum?

Thanks!
Alicia

install.packages(nleqslv,
lib=Users/Alicia/Documents/R.Codes/R.Packages/)

require(nleqslv)

## Global Parameters 
beeta=0.8
pq=1
L=12600
theta=0.6
psale=0.6
mu=psale*(1-theta)
alphah=0.15
Cg=6240
Cs=2820
A= 100
D=0.0001
greekp=0.43
K=10

# Species Parameters ##

b1=0.38
p1=16654
v1 = 0.28
N1=6000
g1=1
delta1=1

b2=0.4
p2=2797
v2 = 0.31
N2=1
g2=1
delta2=1

### Define functions with vector x = c(Lg, Ls, gamma1, gamma2, lamda)

firstordercond - function (x) {

y=numeric(4)


y[1]=(alphah/x[3])-(x[5]*((p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])
+
  delta1*theta*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])))


y[2]=(alphah/x[4])-(x[5]*((p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])
+
  delta2*theta*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])))


y[3]=((alphah*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
+
((alphah*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2]))
+
  x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
(b1*(1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K*((N1/A)*g1^(greekp))*x[1]^(b1-1))
+

(b2*(1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K*((N2/A)*g2^(greekp))*x[1]^(b2-1))
-
(delta1*theta*x[3]*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)) -
(delta2*theta*x[4]*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))-Cg*(((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)+((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))


y[4]=((alphah*(2*v1*N1*D))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
+ ((alphah*(2*v2*N2*D))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])) +
  x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
((1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K*(2*v1*N1*D))
+

((1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K*(2*v2*N2*D))
- (delta1*theta*x[3]*(2*v1*N1*D)) -
  (delta2*theta*x[4]*(2*v2*N2*D)))-Cs*((2*v1*N1*D)+(2*v2*N2*D)))

  result=(x)

}

Xstart=c(1, 200, 0.5, 0.5, 12)
fstart= firstordercond(Xstart)

nleqslv(Xstart, firstordercond)



-- 
Alicia Ellis
Postdoc
Gund Institute for Ecological Economics
University of Vermont
617 Main Street
Burlington, VT  05405
(802) 656-1046
http://www.uvm.edu/~aellis5/
http://entomology.ucdavis.edu/faculty/scott/aellis/

[[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] lme: subject-specific slopes.

2012-12-04 Thread Brian S Cade
John:  I think you want the output from coef(fitRIRT).   The 
ranef(fitRIRT) will give you the subject specific random effect deviations 
from the fixed effects means.  The coef(fitRIRT) will give you the 
combination of the fixed effect means with the subject specific random 
effect deviations and, thus, in the scale you are expecting.

Brian

Brian S. Cade, PhD

U. S. Geological Survey
Fort Collins Science Center
2150 Centre Ave., Bldg. C
Fort Collins, CO  80526-8818

email:  brian_c...@usgs.gov
tel:  970 226-9326



From:
John Sorkin jsor...@grecc.umaryland.edu
To:
r-help@r-project.org, Kenneth Frost kfr...@wisc.edu
Date:
12/04/2012 09:10 AM
Subject:
Re: [R] lme: subject-specific slopes.
Sent by:
r-help-boun...@r-project.org



Ken,
Thank you for your help. ranef(fitRIRT) does not give me what I expect. 
The subject-specific slopes, and subject-specific intercepts are not 
anywhere close to what I would expect them to be; the mean of the 
subject-specfic values should be close to those reported by 
summary(fitRIRT) and they are not. As you will see by examining the 
material below, the subject-specific slopes are off by many order of 
magnitude. The intercepts are also far from the value reported in 
summary(fitRIRT). Do you have any thoughts?
Thanks,
John
 
 fitRIRT - lme(echogen~time,random=~ 
1+time|subject,data=repeatdata,na.action=na.omit) summary(fitRIRT)Linear 
mixed-effects model fit by REML
 Data: repeatdata 
  AIC  BIClogLik
  495.097 507.2491 -241.5485

Random effects:
 Formula: ~1 + time | subject
 Structure: General positive-definite, Log-Cholesky parametrization
StdDev   Corr 
(Intercept) 1.917511e+01 (Intr)
time2.032276e-04 0 
Residual1.044601e+01 

Fixed effects: echogen ~ time 
   Value Std.Error DF   t-value p-value
(Intercept) 64.54864  4.258235 32 15.158543  0.
time 0.35795  0.227080 32  1.576307  0.1248
 Correlation: 
 (Intr)
time -0.242

Standardized Within-Group Residuals:
Min  Q1 Med  Q3 Max 
-1.61362755 -0.52710871  0.02948008  0.41793322  1.77340082 

Number of Observations: 58
Number of Groups: 25  ranef(fitRIRT)   (Intercept)  time
1-3.278112  2.221016e-09
2   -35.400618  4.314995e-08
311.493110 -6.797543e-09
4   -16.209586 -7.070834e-08
5 3.585227 -2.389705e-08
6 1.614320 -1.967700e-09
7 8.346905  5.827094e-08
830.917812 -3.768584e-08
9-0.394101 -9.158251e-09
104.437509 -4.057971e-08
11   31.956597 -2.126275e-08
12   41.567402 -4.853942e-08
13  -10.723993  1.307152e-08
14   -4.554837  5.551908e-09
15   -4.554501  4.815086e-08
16   13.296985 -3.743967e-08
17   -8.255439  1.733238e-08
18  -21.317239  2.203885e-08
19  -13.480159  2.194016e-08
20  -13.044766  2.269168e-08
21   11.639198 -1.418706e-08
22  -27.457388 -1.154099e-08
232.194001 -5.509119e-09
24   -3.992646  7.682188e-08
251.614320 -1.967700e-09
 
 


 
 
John David Sorkin M.D., Ph.D.
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing) 
Kenneth Frost kfr...@wisc.edu 12/4/2012 10:44 AM 
I think this might be close to what you want. ranef(fitRIRT)

On 12/04/12, John Sorkin  wrote:
 I am running a random intercept random slope regression:
 
 fitRIRT - lme(echogen~time,random=~ 
1+time|subject,data=repeatdata,na.action=na.omit)
 summary(fitRIRT)
 
 I would like to get the subject-specific slopes, i.e. the slope that the 
model computes for each subject. If I have 10-subjects I should have 
10-slopes. I don't see the slope when I look at the items listed in 
 names(summary(fitRIRT)
 nor when I look at the items listed in 
 names(fitRIRT)
 
 Thanks 
 John
 
 
 
 
 
 
 John David Sorkin M.D., Ph.D.
 Chief, Biostatistics and Informatics
 University of Maryland School of Medicine Division of Gerontology
 Baltimore VA Medical Center
 10 North Greene Street
 GRECC (BT/18/GR)
 Baltimore, MD 21201-1524
 (Phone) 410-605-7119
 (Fax) 410-605-7913 (Please call phone number above prior to faxing)
 Confidentiality Statement:
 This email message, including any attachments, is for the sole use of 
the intended recipient(s) and may contain confidential and privileged 
information. Any unauthorized use, disclosure or distribution is 
prohibited. If you are not the intended recipient, please contact the 
sender by reply email and destroy all copies of the original message. 
 
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

--
Kenneth Frost
Graduate Research Assistant - Dept. of Plant Pathology
University of Wisconsin - Madison
Lab: 

Re: [R] computing marginal values based on multiple columns?

2012-12-04 Thread arun
HI,

I am not sure the output you wanted is correct: 


sample1 sample2 sample3
1      1.0     0     0.5


because
0.2*colMeans(x[,-4])
sample1 sample2 sample3 
#  28.40   24.08   21.36 


This might help you:
apply(x[-4],2,function(y) length(y[y 0.2*mean(y)  
x$class==a])/length(x[x$class==a]))
#sample1 sample2 sample3 
  #  0.0 0.0 0.5 
A.K.



- Original Message -
From: Simon simonzm...@gmail.com
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 4:49 AM
Subject: [R] computing marginal values based on multiple columns?

Hello all,

I have what feels like a simple problem, but I can't find an simple
answer. Consider this data frame:

 x - data.frame(sample1=c(35,176,182,193,124),
sample2=c(198,176,190,23,15), sample3=c(12,154,21,191,156),
class=c('a','a','c','b','c'))

 x
  sample1 sample2 sample3 class
1      35     198      12     a
2     176     176     154     a
3     182     190      21     c
4     193      23     191     b
5     124      15     156     c

Now I wish to know: for each sample, for values  20% of the sample mean,
what percentage of those are class a?

I want to end up with a table like:

   sample1 sample2 sample3
1      1.0     0     0.5

I can calculate this for an individual sample using this rather clumsy
expression:

length(which(x$sample1  mean(x$sample1)  x$class=='a')) /
length(which(x$sample1  mean(x$sample1)))

I'd normally propagate it across the data frame using apply, but I
can't because it depends on more than one column.

Any help much appreciated!

Cheers,

Simon

    [[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] Help for a function

2012-12-04 Thread Jeremy Miles
What are you expecting?

What do you get?

What is the problem?

J

On 4 December 2012 06:01, anoumou teko_maur...@yahoo.fr wrote:
 Hello all,
 I need a help.
 I am modeling a disease and a create a R function like that:

 Lambda-function (x,date1,r,h,a){
   ndate1 - as.Date(date1, %d/%m/%Y)
   t1 - as.numeric(ndate1)
   x[order(x$i),]
   t -x[,t]
   i -x[,i]
   CONTAGIEUX -x[,CONTAGIEUX]
   while ( t1  min(t) ){
   for (i in 1:length(i) ){
 {for (j in 1:CONTAGIEUX[length(CONTAGIEUX)]){
   res1[j] -(a*h)
   res2 -sum( res1[j])
 }
  }
   lambda[i] - r*res2
   }
   }
 x-data.frame(x,lambda)
 x
 }

 on such data :

 DATEi   Symptomes   t   Incubation  CONTAGIEUX
 1   2009-04-29   Canada 13  14363   13  13
 2   2009-05-01   Israel 2   14365   2   2
 3   2009-05-09  argentina   1   14373   1   1
 5   2009-05-09  australia   1   14373   1   1
 6   2009-05-10  australia   1   14374   2   2
 7   2009-04-29  Austria 1   14363   1   1
 8   2009-04-30  Austria 1   14364   2   2
 9   2009-05-01  Austria 1   14365   2   3
 10  2009-05-02  Austria 1   14366   2   4
 11  2009-05-03  Austria 1   14367   2   5
 17  2009-05-09  Austria 1   14373   2   7
 18  2009-05-10  Austria 1   14374   2   7
 19  2009-05-08  brasil  4   14372   4   4
 20  2009-05-09  brazil  6   14373   6   6
 21  2009-05-10  brazil  6   14374   12  12
 22  2009-05-02  canada  51  14366   51  51
 23  2009-05-03  canada  85  14367   136 136
 24  2009-05-04  canada  101 14368   186 237
 31  2009-04-27  Canada  6   14361   6   6
 32  2009-04-28  Canada  6   14362   6   6
 33  2009-04-30  Canada  19  14364   25  25
 34  2009-05-01  Canada  34  14365   53  59
 35  2009-05-01  China,HongKong, SAR 1   14365   1   1
 36  2009-05-02  China,HongKong, SAR 1   14366   2   2
 37  2009-05-03  China,HongKong, SAR 1   14367   2   3
 38  2009-05-04  China,HongKong, SAR 1   14368   2   4
 44  2009-05-10  China,HongKong, SAR 1   14374   2   7
 45  2009-05-04  colombia1   14368   1   1
 46  2009-05-05  colombia1   14369   2   2
 47  2009-05-06  colombia

 But i do not get the results,i try by all means but i d'ont understant the
 problem.
 Thanks for your help.



 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Help-for-a-function-tp4652054.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] Labelling x axis in plot function

2012-12-04 Thread T Bal
Hi,
In the plot function I want to label x axis as the numbers between 1 and 12
(so 1, 2, 3, 4, 5, ..., 12). How should I do it? The range of x values are
different than this range. Thanks!


Kind regards,
T. Bal

[[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] Labelling x axis in plot function

2012-12-04 Thread Rui Barradas

Hello,

You should provide sample data and code.

plot(3:10, xlim = c(0, 12), xaxt = n)
axis(1, at = 1:12)


See the help page ?par for a description of the graphical parameters 
'xlim' and 'xaxt', and of ?axis.


Hope this helps,

Rui Barradas
Em 04-12-2012 17:24, T Bal escreveu:

Hi,
In the plot function I want to label x axis as the numbers between 1 and 12
(so 1, 2, 3, 4, 5, ..., 12). How should I do it? The range of x values are
different than this range. Thanks!


Kind regards,
T. Bal

[[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] Solve system of equations (nleqslv) only returns origin

2012-12-04 Thread Berend Hasselman

On 04-12-2012, at 17:06, Alicia Ellis wrote:

 I'm solving 4 complex equations simultaneously.  Code is below.  The code
 returns only zero's for the solution though there should also be a non-zero
 result.  I'm pretty confident that the equations are correct because they
 are straight from a published paper and I checked them pretty thoroughly.
 The parameter values I used are from the published paper as well.  Any
 suggestions for how I can find the maximum non-zero solution instead of
 going straight to the minimum?
 

Are you trying to minimize a function by solving the first order condition?
If yes then you should try functions such optim, nlminb and there are many more.

See below for more comments.

 Thanks!
 Alicia
 
 install.packages(nleqslv,
 lib=Users/Alicia/Documents/R.Codes/R.Packages/)
 
 require(nleqslv)
 
 ## Global Parameters 
 beeta=0.8
 pq=1
 L=12600
 theta=0.6
 psale=0.6
 mu=psale*(1-theta)
 alphah=0.15
 Cg=6240
 Cs=2820
 A= 100
 D=0.0001
 greekp=0.43
 K=10
 
 # Species Parameters ##
 
 b1=0.38
 p1=16654
 v1 = 0.28
 N1=6000
 g1=1
 delta1=1
 
 b2=0.4
 p2=2797
 v2 = 0.31
 N2=1
 g2=1
 delta2=1
 
 ### Define functions with vector x = c(Lg, Ls, gamma1, gamma2, lamda)
 
 firstordercond - function (x) {
 
y=numeric(4)
 
 
 y[1]=(alphah/x[3])-(x[5]*((p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])
 +
  delta1*theta*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])))
 
 
 y[2]=(alphah/x[4])-(x[5]*((p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])
 +
  delta2*theta*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])))
 
 
 y[3]=((alphah*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
 +
 ((alphah*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2]))
 +
  x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
 (b1*(1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K*((N1/A)*g1^(greekp))*x[1]^(b1-1))
 +
 
 (b2*(1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K*((N2/A)*g2^(greekp))*x[1]^(b2-1))
 -
(delta1*theta*x[3]*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)) -
 (delta2*theta*x[4]*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))-Cg*(((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)+((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))
 
 
 y[4]=((alphah*(2*v1*N1*D))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
 + ((alphah*(2*v2*N2*D))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])) +
  x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
 ((1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K*(2*v1*N1*D))
 +
 
 ((1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K*(2*v2*N2*D))
 - (delta1*theta*x[3]*(2*v1*N1*D)) -
  (delta2*theta*x[4]*(2*v2*N2*D)))-Cs*((2*v1*N1*D)+(2*v2*N2*D)))
 
  result=(x)
 

what is this statement supposed to do?
You are actually returning the input.
And why like that?
Just x or return(x) is quite sufficient.

You should return the vector y i.e. function values.
But y has length 4 and x has length 4.

So where is the fifth value for y?

}
 
 Xstart=c(1, 200, 0.5, 0.5, 12)
 fstart= firstordercond(Xstart)
 

If you print fstart you will see that it is identical to Xstart.

You need to rethink you firstordercond() function.
It's not correct.

Berend

 nleqslv(Xstart, firstordercond)
 
 
 
 -- 
 Alicia Ellis
 Postdoc
 Gund Institute for Ecological Economics
 University of Vermont
 617 Main Street
 Burlington, VT  05405
 (802) 656-1046
 http://www.uvm.edu/~aellis5/
 http://entomology.ucdavis.edu/faculty/scott/aellis/
 
   [[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] odd behavior of browser()

2012-12-04 Thread David Romano
Hi everyone,

I normally include a call to browser() as I'm working out the kinks in my
scripts, and I am always able to step through each line by hitting
Return, but for some reason, in the scripts I'm working on now, hitting
Return seems to cause execution of *all* the lines in my script.  I've
restarted R several times in case it was stuck in a bad state for some
reason, but I'm consistently getting this behavior anyway.  Has anyone run
into this problem before?  Maybe I inadvertently reset preferences?

An example which produces this behavior is the following:

file bugcheck.r:

browser()

a - 1
b - 2

 source(bugcheck.r)
Called from: eval(expr, envir, enclos)
Browse[1]  Return

 ls()
[1] a b
 a
[1] 1
 b
[1] 2

I'd be grateful for any help in resolving this!

Thanks,
David Romano

[[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] Solve system of equations (nleqslv) only returns origin

2012-12-04 Thread Berend Hasselman

On 04-12-2012, at 18:50, Berend Hasselman wrote:

 
 
 You should return the vector y i.e. function values.
 But y has length 4 and x has length 4.
 

x has length 5 of course.

Berend

 So where is the fifth value for y?
 
   }
 
 Xstart=c(1, 200, 0.5, 0.5, 12)
 fstart= firstordercond(Xstart)
 
 
 If you print fstart you will see that it is identical to Xstart.
 
 You need to rethink you firstordercond() function.
 It's not correct.
 
 Berend
 
 nleqslv(Xstart, firstordercond)
 
 
 
 -- 
 Alicia Ellis
 Postdoc
 Gund Institute for Ecological Economics
 University of Vermont
 617 Main Street
 Burlington, VT  05405
 (802) 656-1046
 http://www.uvm.edu/~aellis5/
 http://entomology.ucdavis.edu/faculty/scott/aellis/
 
  [[alternative HTML version deleted]]
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

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


[R] error reading xlsm file with read.xls

2012-12-04 Thread Sebastian Kruk
Dear all,

I cannot reading a .xlsm file using read.xls.

I executed:

read.xls(resultados.xlsm,
colNames = TRUE,
sheet = 1,
type = data.frame,
from = 1,
rowNames = NA,
colClasses = character,
checkNames = TRUE,
dateTime = numeric,
naStrings = NA,
stringsAsFactors = F)
Error:

Call(ReadXls, file, colNames, sheet, type, from, rowNames, :
Incorrect number of arguments (11), expecting 10 for 'ReadXls'
If I just write

read.xls(resultados.xlsm)
It give me the same error.

Regards,

Sebastián.

[[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] odd behavior of browser()

2012-12-04 Thread R. Michael Weylandt
Untested, I think it's the blank line in your script which exits the
debugger and then you're seeing regular execution.

MW

On Tue, Dec 4, 2012 at 5:54 PM, David Romano drom...@stanford.edu wrote:
 Hi everyone,

 I normally include a call to browser() as I'm working out the kinks in my
 scripts, and I am always able to step through each line by hitting
 Return, but for some reason, in the scripts I'm working on now, hitting
 Return seems to cause execution of *all* the lines in my script.  I've
 restarted R several times in case it was stuck in a bad state for some
 reason, but I'm consistently getting this behavior anyway.  Has anyone run
 into this problem before?  Maybe I inadvertently reset preferences?

 An example which produces this behavior is the following:

 file bugcheck.r:

 browser()

 a - 1
 b - 2

 source(bugcheck.r)
 Called from: eval(expr, envir, enclos)
 Browse[1]  Return

 ls()
 [1] a b
 a
 [1] 1
 b
 [1] 2

 I'd be grateful for any help in resolving this!

 Thanks,
 David Romano

 [[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] Read.csv

2012-12-04 Thread Torus Insurance
Hi list,
I am using read.csv to read data from csf files, but noticed that the
numeric data (those larger than 10 power 9) are rounded to the nearest
million (10 power 6). Any solution?
Thanks
Arvin

-- 
Sent from my mobile device

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

2012-12-04 Thread Duncan Murdoch

On 04/12/2012 1:02 PM, Torus Insurance wrote:

Hi list,
I am using read.csv to read data from csf files, but noticed that the
numeric data (those larger than 10 power 9) are rounded to the nearest
million (10 power 6). Any solution?


What makes you think that is happening?  R rounds values for printing, 
but generally maintains about 15 digits internally.


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] Read.csv

2012-12-04 Thread Sarah Goslee
Hi Arvin,

How are you preparing the data to be read: a spreadsheet?
How are you reading the data?

How have you verified that the CSV file is correct?
How have you verified that the data frame is incorrect?

Can you provide a reproducible example using a small portion of your dataset?

Sarah

On Tue, Dec 4, 2012 at 1:02 PM, Torus Insurance
torusinsuran...@gmail.com wrote:
 Hi list,
 I am using read.csv to read data from csf files, but noticed that the
 numeric data (those larger than 10 power 9) are rounded to the nearest
 million (10 power 6). Any solution?
 Thanks
 Arvin

 --
 Sent from my mobile device


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

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


[R] monte carlo simulation on R

2012-12-04 Thread Adel ESSAFI
Hello,

How can I make a monte carlo simulation on R?

Regards
Adel


-- 
PhD candidate in Computer Science
Address
3 avenue lamine, cité ezzahra, Sousse 4000
Tunisia
tel: +216 97 246 706 (+33640302046 jusqu'au 15/6)
fax: +216 71 391 166

[[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] odd behavior of browser()

2012-12-04 Thread Duncan Murdoch

On 04/12/2012 12:54 PM, David Romano wrote:

Hi everyone,

I normally include a call to browser() as I'm working out the kinks in my
scripts, and I am always able to step through each line by hitting
Return, but for some reason, in the scripts I'm working on now, hitting
Return seems to cause execution of *all* the lines in my script.  I've
restarted R several times in case it was stuck in a bad state for some
reason, but I'm consistently getting this behavior anyway.  Has anyone run
into this problem before?  Maybe I inadvertently reset preferences?


I wouldn't have expected that to work.  Calling browser() from within a 
function will let you step through the function, but calling it from 
within a script doesn't.  Do you really have some scripts where this worked?


Duncan Murdoch




An example which produces this behavior is the following:

file bugcheck.r:

browser()

a - 1
b - 2

 source(bugcheck.r)
Called from: eval(expr, envir, enclos)
Browse[1]  Return

 ls()
[1] a b
 a
[1] 1
 b
[1] 2

I'd be grateful for any help in resolving this!

Thanks,
David Romano

[[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] monte carlo simulation on R

2012-12-04 Thread R. Michael Weylandt
replicate(1000, sum(rnorm(50)^2-rchisq(50, 3)))

Or you know, many other things...

Michael

On Tuesday, December 4, 2012, Adel ESSAFI wrote:

 Hello,

 How can I make a monte carlo simulation on R?

 Regards
 Adel


 --
 PhD candidate in Computer Science
 Address
 3 avenue lamine, cité ezzahra, Sousse 4000
 Tunisia
 tel: +216 97 246 706 (+33640302046 jusqu'au 15/6)
 fax: +216 71 391 166

 [[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] Winbugs from R

2012-12-04 Thread Bob O'Hara
This is a BUGS problem, not an R problem, so would be better asked on the
BUGS list (a practical hint: if it's not working through R, run the code
directly in BUGS: the error reporting is better). Anyway, here BUGS manages
to tell R to tell you what's wrong:

variable n is not defined

And looking at your data, it's right. n is not defined


On 4 December 2012 14:17, Günal Bilek gunal...@hotmail.com wrote:

 Hi,
 I am trying to covert a Winbugs code into R code. Here is the winbugs code
 model{# model’s likelihoodfor (i in 1:n){time[i] ~ dnorm( mu[i], tau ) #
 stochastic componenent# link and linear predictormu[i] - beta0 + beta1 *
 cases[i] + beta2 * distance[i]}# prior distributionstau ~ dgamma( 0.01,
 0.01 )beta0 ~ dnorm( 0.0, 1.0E-4)beta1 ~ dnorm( 0.0, 1.0E-4)beta2 ~ dnorm(
 0.0, 1.0E-4)# definition of sigmas2-1/taus -sqrt(s2)# calculation of the
 sample variancefor (i in 1:n){ c.time[i]-time[i]-mean(time[]) }sy2 -
 inprod( c.time[], c.time[] )/(n-1)# calculation of Bayesian version R
 squaredR2B - 1 - s2/sy2# Expected y for a typical delivery timetypical.y
 - beta0 + beta1 * mean(cases[]) + beta2 * mean(distance[])}INITSlist(
 tau=1, beta0=1, beta1=0, beta2=0 )DATA (LIST)list( n=25,time = c(16.68,
 11.5, 12.03, 14.88, 13.75, 18.11, 8, 17.83,79.24, 21.5, 40.33, 21, 13.5,
 19.75, 24, 29, 15.35,19, 9.5, 35.1, 17.9, 52.32, 18.75, 19.83,
 10.75),distance = c(560, 220, 340, 80, 150, 330, 110, 210, 1460,605, 688,
 215, 255, 462, 448, 776, 200, 132,36, 770, 140, 810, 450, 635, 150),cases =
 c( 7, 3, 3, 4, 6, 7, 2, 7, 30, 5, 16, 10, 4, 6, 9,10, 6, 7, 3, 17, 10, 26,
 9, 8, 4) )

 I want to do this in R. So, I copied the model and pasted into a txt file
 named reg. Here is the R code I used
 time -
 c(16.68,11.5,12.03,14.88,13.75,18.11,8,17.83,79.24,21.5,40.33,21,13.5,19.75,24,29,15.35,19,9.5,35.1,17.9,52.32,18.75,19.83,10.75)cases
 - c(7,3,3,4,6,7,2,7,30,5,16,10,4,6,9,10,6,7,3,17,10,26,9,8,4)distance -
 c(560,220,340,80,150,330,110,210,1260,605,688,215,255,462,448,776,200,132,36,770,140,810,450,635,150)
 data - list(time,cases,distance)
 inits - function(){list(tau=1,beta1=0,beta2=0,beta3=0)}
 sim - bugs(data, inits, model.file =
 C:/Users/Gunal/Desktop/dummy/reg.txt,parameters = c(beta1,
 beta2,beta3),n.chains = 3, n.iter = 1000,bugs.directory =
 D:/PROGRAMLAR/WinBUGS14/,debug=TRUE)

 Winbugs is producing this error page
 display(log)check(C:/Users/Gunal/Desktop/dummy/reg.txt)model is
 syntactically
 correctdata(C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/data.txt)data
 loadedcompile(3)variable n is not
 definedinits(1,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/inits1.txt)command
 #Bugs:inits cannot be executed (is greyed
 out)inits(2,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/inits2.txt)command
 #Bugs:inits cannot be executed (is greyed
 out)inits(3,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/inits3.txt)command
 #Bugs:inits cannot be executed (is greyed out)gen.inits()command
 #Bugs:gen.inits cannot be executed (is greyed
 out)thin.updater(1)update(500)command #Bugs:update cannot be executed (is
 greyed out)set(beta1)command #Bugs:set cannot be executed (is greyed
 out)set(beta2)command #Bugs:set cannot be executed (is greyed
 out)set(beta3)command #Bugs:set cannot be executed (is greyed
 out)set(deviance)command #Bugs:set cannot be executed (is greyed
 out)dic.set()command #Bugs:dic.set cannot be executed (is greyed
 out)update(500)command #Bugs:update cannot be executed (is greyed
 out)coda(*,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/coda)command
 #Bugs:coda cannot be executed (is greyed out)stats(*)command #Bugs:stats
 cannot be executed (is greyed out)dic.stats()
 DIChistory(*,C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/history.odc)command
 #Bugs:history cannot be executed (is greyed
 out)save(C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/log.odc)save(C:/Users/Gunal/AppData/Local/Temp/RtmpCYxtDZ/log.txt)

 Any help would be greatly appreciated.
 Cheers
 Gunal
 [[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.




-- 
Bob O'Hara

Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany

Tel: +49 69 798 40226
Mobile: +49 1515 888 5440
WWW:   http://www.bik-f.de/root/index.php?page_id=219
Blog: http://occamstypewriter.org/boboh/
Journal of Negative Results - EEB: www.jnr-eeb.org

[[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] Large loops in R

2012-12-04 Thread Charles Novaes de Santana
Dear all,

I need to access data from a large matrix (48000 x 48000) and to do it
I am trying to run two loops using for command. Surely it is been a
very slow job.

I heard that for is not the best option to perform large loops in R,
but I don't really know what would be the best (fast) option. sapply?
vapply? Could anyone help me with this issue, please?

Thank you very much for your attention and for any help!

Best regards,

Charles

-- 
Um axé! :)

--
Charles Novaes de Santana
http://www.imedea.uib-csic.es/~charles
PhD student - Global Change
Laboratorio Internacional de Cambio Global
Department of Global Change Research
Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
Calle Miquel Marques 21, 07190
Esporles - Islas Baleares - España

Office phone - +34 971 610 896
Cell phone - +34 660 207 940

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

2012-12-04 Thread R. Michael Weylandt michael.weyla...@gmail.com


On Dec 4, 2012, at 6:47 PM, Charles Novaes de Santana 
charles.sant...@gmail.com wrote:

 Dear all,
 
 I need to access data from a large matrix (48000 x 48000) and to do it
 I am trying to run two loops using for command. Surely it is been a
 very slow job.
 
 I heard that for is not the best option to perform large loops in R,
 but I don't really know what would be the best (fast) option. sapply?
 vapply? Could anyone help me with this issue, please?

What exactly are you trying to do? It's likely doable with a few vectorized 
operations. 

Michael 


 
 Thank you very much for your attention and for any help!
 
 Best regards,
 
 Charles
 
 -- 
 Um axé! :)
 
 --
 Charles Novaes de Santana
 http://www.imedea.uib-csic.es/~charles
 PhD student - Global Change
 Laboratorio Internacional de Cambio Global
 Department of Global Change Research
 Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
 Calle Miquel Marques 21, 07190
 Esporles - Islas Baleares - España
 
 Office phone - +34 971 610 896
 Cell phone - +34 660 207 940
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] control point size of superscript when labeling axes with title()

2012-12-04 Thread Chris Solomon
Hi-

A journal has asked me to make all of my text annotations on a figure at 
10-point size. For the most part this is easy, e.g. by creating figures with:
   pdf(..., family='Times', pointsize=10)

But where I have superscripts (or subscripts) in axis labels, the default seems 
to be to shrink the superscripted text slightly. For example this code:
   title(ylab=expression(paste('Respiration (mg  ',O[2],'  ',L^-1,' 
',d^-1,')',sep=' ')),outer=T,line=0.3)
produces superscripted numbers at approximately 7 point.

I have been poking around for a solution but not having much luck. The 
textstyle() function keeps superscripted text at original size if you use it 
within a text() call, but I can't figure out an equivalent solution within a 
title() call.

Thanks for any suggestions
Chris

--
Dr. Chris Solomon
Assistant Professor
Dept. of Natural Resource Sciences
McGill University




[[alternative HTML version deleted]]

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


Re: [R] Large loops in R

2012-12-04 Thread Charles Novaes de Santana
Dear Michael,

Thank you for your answer.

I have 2 matrices. Each position of the matrices is a weight. And I
need to calculate the following sum of differences:

Considering:
mat1 and mat2 - two matrices (each of them 48000 x 48000).
d1 and d2 - two constant values.

sum-0;
for(i in 1:nrows1){
for(j in 1:nrows2){
sum-sum+ ( ( (mat1(i,j)/d1) -
(mat2(i,j)/d2) )^2 )
}
}
}

I was wondering if there is a better way to do this sum.

Thank you for your attention!

Best,

Charles

On Tue, Dec 4, 2012 at 7:54 PM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:


 On Dec 4, 2012, at 6:47 PM, Charles Novaes de Santana 
 charles.sant...@gmail.com wrote:

 Dear all,

 I need to access data from a large matrix (48000 x 48000) and to do it
 I am trying to run two loops using for command. Surely it is been a
 very slow job.

 I heard that for is not the best option to perform large loops in R,
 but I don't really know what would be the best (fast) option. sapply?
 vapply? Could anyone help me with this issue, please?

 What exactly are you trying to do? It's likely doable with a few vectorized 
 operations.

 Michael



 Thank you very much for your attention and for any help!

 Best regards,

 Charles

 --
 Um axé! :)

 --
 Charles Novaes de Santana
 http://www.imedea.uib-csic.es/~charles
 PhD student - Global Change
 Laboratorio Internacional de Cambio Global
 Department of Global Change Research
 Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
 Calle Miquel Marques 21, 07190
 Esporles - Islas Baleares - España

 Office phone - +34 971 610 896
 Cell phone - +34 660 207 940

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



-- 
Um axé! :)

--
Charles Novaes de Santana
http://www.imedea.uib-csic.es/~charles
PhD student - Global Change
Laboratorio Internacional de Cambio Global
Department of Global Change Research
Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
Calle Miquel Marques 21, 07190
Esporles - Islas Baleares - España

Office phone - +34 971 610 896
Cell phone - +34 660 207 940

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

2012-12-04 Thread Sarah Goslee
Without a reproducible example it's hard to tell for certain, but what
about simply (assuming nrows2 is actually columns):

sum((mat1/d1 - mat2/d2)^2)

R is smart enough to understand elementwise manipulation of a matrix:
you shouldn't need a loop at all.

Sarah

On Tue, Dec 4, 2012 at 2:27 PM, Charles Novaes de Santana
charles.sant...@gmail.com wrote:
 Dear Michael,

 Thank you for your answer.

 I have 2 matrices. Each position of the matrices is a weight. And I
 need to calculate the following sum of differences:

 Considering:
 mat1 and mat2 - two matrices (each of them 48000 x 48000).
 d1 and d2 - two constant values.

 sum-0;
 for(i in 1:nrows1){
 for(j in 1:nrows2){
 sum-sum+ ( ( (mat1(i,j)/d1) -
 (mat2(i,j)/d2) )^2 )
 }
 }
 }

 I was wondering if there is a better way to do this sum.

 Thank you for your attention!

 Best,

 Charles

 On Tue, Dec 4, 2012 at 7:54 PM, R. Michael Weylandt
 michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:


 On Dec 4, 2012, at 6:47 PM, Charles Novaes de Santana 
 charles.sant...@gmail.com wrote:

 Dear all,

 I need to access data from a large matrix (48000 x 48000) and to do it
 I am trying to run two loops using for command. Surely it is been a
 very slow job.

 I heard that for is not the best option to perform large loops in R,
 but I don't really know what would be the best (fast) option. sapply?
 vapply? Could anyone help me with this issue, please?

 What exactly are you trying to do? It's likely doable with a few vectorized 
 operations.

 Michael


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

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


Re: [R] Large loops in R

2012-12-04 Thread Peter Langfelder
On Tue, Dec 4, 2012 at 11:27 AM, Charles Novaes de Santana
charles.sant...@gmail.com wrote:
 Dear Michael,

 Thank you for your answer.

 I have 2 matrices. Each position of the matrices is a weight. And I
 need to calculate the following sum of differences:

 Considering:
 mat1 and mat2 - two matrices (each of them 48000 x 48000).
 d1 and d2 - two constant values.

 sum-0;
 for(i in 1:nrows1){
 for(j in 1:nrows2){
 sum-sum+ ( ( (mat1(i,j)/d1) -
 (mat2(i,j)/d2) )^2 )
 }
 }
 }

 I was wondering if there is a better way to do this sum.

sum( (mat1/d1-mat2/d2)^2)

Correct me if I'm wrong though - aren't matrices of 48x times 48k
larger than what R can handle at present?

HTH

Peter

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

2012-12-04 Thread Sam Steingold
How do I convert a list to a matrix?

--8---cut here---start-8---
list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17), 
c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25), 
c(45, 19), c(5e+05, 16))
as.matrix(a)
  [,1] 
 [1,] Numeric,2
 [2,] Numeric,2
 [3,] Numeric,2
 [4,] Numeric,2
 [5,] Numeric,2
 [6,] Numeric,2
 [7,] Numeric,2
 [8,] Numeric,2
 [9,] Numeric,2
--8---cut here---end---8---

thanks!

-- 
Sam Steingold (http://sds.podval.org/) on Ubuntu 12.04 (precise) X 11.0.11103000
http://www.childpsy.net/ http://palestinefacts.org http://dhimmi.com
http://jihadwatch.org http://www.PetitionOnline.com/tap12009/ http://memri.org
Rhinoceros has poor vision, but, due to his size, it's not his problem.

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

2012-12-04 Thread Suzen, Mehmet
There is a book on MC  with R:

Introducing Monte Carlo Methods with R by Robert/Casella:
http://www.springer.com/statistics/computational+statistics/book/978-1-4419-1575-7

On 4 December 2012 19:21, Adel ESSAFI adeless...@gmail.com wrote:
 Hello,

 How can I make a monte carlo simulation on R?

 Regards
 Adel


 --
 PhD candidate in Computer Science
 Address
 3 avenue lamine, cité ezzahra, Sousse 4000
 Tunisia
 tel: +216 97 246 706 (+33640302046 jusqu'au 15/6)
 fax: +216 71 391 166

 [[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] Large loops in R

2012-12-04 Thread Charles Novaes de Santana
On Tue, Dec 4, 2012 at 8:43 PM, Peter Langfelder
peter.langfel...@gmail.com wrote:
 On Tue, Dec 4, 2012 at 11:27 AM, Charles Novaes de Santana
 charles.sant...@gmail.com wrote:
 Dear Michael,

 Thank you for your answer.

 I have 2 matrices. Each position of the matrices is a weight. And I
 need to calculate the following sum of differences:

 Considering:
 mat1 and mat2 - two matrices (each of them 48000 x 48000).
 d1 and d2 - two constant values.

 sum-0;
 for(i in 1:nrows1){
 for(j in 1:nrows2){
 sum-sum+ ( ( (mat1(i,j)/d1) -
 (mat2(i,j)/d2) )^2 )
 }
 }
 }

 I was wondering if there is a better way to do this sum.

 sum( (mat1/d1-mat2/d2)^2)

 Correct me if I'm wrong though - aren't matrices of 48x times 48k
 larger than what R can handle at present?

 HTH

 Peter

hmmm I didn't know that the limitation of R was below this value. I
found this error message:

Error in matrix(0, 48000, 48000) : too many elements specified

but I thought it was a machine limitation (and I was asking for access
to a better machine in my labs...). Thanks for clarifying it.

Well, when Sarah gave me the answer for my problem, I got a new one :)
Thank you, Sarah and Peter.

So, is there any other way to trick R and allocate such large matrices?

Best,

Charles
-- 
Um axé! :)

--
Charles Novaes de Santana
http://www.imedea.uib-csic.es/~charles
PhD student - Global Change
Laboratorio Internacional de Cambio Global
Department of Global Change Research
Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
Calle Miquel Marques 21, 07190
Esporles - Islas Baleares - España

Office phone - +34 971 610 896
Cell phone - +34 660 207 940

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

2012-12-04 Thread Mark Lamias
Try:

matrix(unlist(a), ncol=2, byrow=T)

--Mark Lamias





 From: Sam Steingold s...@gnu.org
To: r-help@r-project.org 
Sent: Tuesday, December 4, 2012 3:09 PM
Subject: [R] list to matrix?

How do I convert a list to a matrix?

--8---cut here---start-8---
list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17), 
    c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25), 
    c(45, 19), c(5e+05, 16))
as.matrix(a)
      [,1]    
[1,] Numeric,2
[2,] Numeric,2
[3,] Numeric,2
[4,] Numeric,2
[5,] Numeric,2
[6,] Numeric,2
[7,] Numeric,2
[8,] Numeric,2
[9,] Numeric,2
--8---cut here---end---8---

thanks!

-- 
Sam Steingold (http://sds.podval.org/) on Ubuntu 12.04 (precise) X 11.0.11103000
http://www.childpsy.net/ http://palestinefacts.org http://dhimmi.com
http://jihadwatch.org http://www.PetitionOnline.com/tap12009/ http://memri.org
Rhinoceros has poor vision, but, due to his size, it's not his problem.

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

2012-12-04 Thread R. Michael Weylandt
On Tue, Dec 4, 2012 at 8:09 PM, Sam Steingold s...@gnu.org wrote:
 How do I convert a list to a matrix?

 --8---cut here---start-8---
 list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17),
 c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25),
 c(45, 19), c(5e+05, 16))

do.call(rbind, LIST)

or

do.call(cbind, LIST)

depending on your desired direction. The do.call syntax takes a
function name and uses a list as arguments to that function, returning
the result. Super useful for these situations where you collect things
and something like

cbind(x[[1]],x[[2]], x[[3]]...)

isn't feasible or possible.

MW

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

2012-12-04 Thread Charles Novaes de Santana
Thank you, Sarah! It is a wonderful new!!! :)

Now I need to solve the other question hehe How to allocate such large matrix :)

best,

Charles

On Tue, Dec 4, 2012 at 8:39 PM, Sarah Goslee sarah.gos...@gmail.com wrote:
 Without a reproducible example it's hard to tell for certain, but what
 about simply (assuming nrows2 is actually columns):

 sum((mat1/d1 - mat2/d2)^2)

 R is smart enough to understand elementwise manipulation of a matrix:
 you shouldn't need a loop at all.

 Sarah

 On Tue, Dec 4, 2012 at 2:27 PM, Charles Novaes de Santana
 charles.sant...@gmail.com wrote:
 Dear Michael,

 Thank you for your answer.

 I have 2 matrices. Each position of the matrices is a weight. And I
 need to calculate the following sum of differences:

 Considering:
 mat1 and mat2 - two matrices (each of them 48000 x 48000).
 d1 and d2 - two constant values.

 sum-0;
 for(i in 1:nrows1){
 for(j in 1:nrows2){
 sum-sum+ ( ( (mat1(i,j)/d1) -
 (mat2(i,j)/d2) )^2 )
 }
 }
 }

 I was wondering if there is a better way to do this sum.

 Thank you for your attention!

 Best,

 Charles

 On Tue, Dec 4, 2012 at 7:54 PM, R. Michael Weylandt
 michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:


 On Dec 4, 2012, at 6:47 PM, Charles Novaes de Santana 
 charles.sant...@gmail.com wrote:

 Dear all,

 I need to access data from a large matrix (48000 x 48000) and to do it
 I am trying to run two loops using for command. Surely it is been a
 very slow job.

 I heard that for is not the best option to perform large loops in R,
 but I don't really know what would be the best (fast) option. sapply?
 vapply? Could anyone help me with this issue, please?

 What exactly are you trying to do? It's likely doable with a few vectorized 
 operations.

 Michael


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



-- 
Um axé! :)

--
Charles Novaes de Santana
http://www.imedea.uib-csic.es/~charles
PhD student - Global Change
Laboratorio Internacional de Cambio Global
Department of Global Change Research
Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
Calle Miquel Marques 21, 07190
Esporles - Islas Baleares - España

Office phone - +34 971 610 896
Cell phone - +34 660 207 940

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

2012-12-04 Thread R. Michael Weylandt
On Tue, Dec 4, 2012 at 8:14 PM, Charles Novaes de Santana
charles.sant...@gmail.com wrote:
 Error in matrix(0, 48000, 48000) : too many elements specified

 but I thought it was a machine limitation (and I was asking for access
 to a better machine in my labs...). Thanks for clarifying it.

 Well, when Sarah gave me the answer for my problem, I got a new one :)
 Thank you, Sarah and Peter.

 So, is there any other way to trick R and allocate such large matrices?


Either

1) Use the development version of R which has large-vector support
(matrices are vectors under the hood)

or

2) A large matrix is usually sparse in structure: use one of the many
sparse representation packages (e.g., Matrix) available.

MW

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

2012-12-04 Thread Sarah Goslee
I don't think there's any reason for the calculation you're doing that
you must have the whole matrix in memory, is there?

Unless there's something more than what you've shown us, you're just
taking the sum of elementwise operations. You can read the matrix in
in manageable chunks, take the sum of that chunk and save the single
value. Repeat, then add them all up at the end.

Sarah

On Tue, Dec 4, 2012 at 3:14 PM, Charles Novaes de Santana
charles.sant...@gmail.com wrote:
 On Tue, Dec 4, 2012 at 8:43 PM, Peter Langfelder
 peter.langfel...@gmail.com wrote:
 On Tue, Dec 4, 2012 at 11:27 AM, Charles Novaes de Santana
 charles.sant...@gmail.com wrote:
 Dear Michael,

 Thank you for your answer.

 I have 2 matrices. Each position of the matrices is a weight. And I
 need to calculate the following sum of differences:

 Considering:
 mat1 and mat2 - two matrices (each of them 48000 x 48000).
 d1 and d2 - two constant values.

 sum-0;
 for(i in 1:nrows1){
 for(j in 1:nrows2){
 sum-sum+ ( ( (mat1(i,j)/d1) -
 (mat2(i,j)/d2) )^2 )
 }
 }
 }

 I was wondering if there is a better way to do this sum.

 sum( (mat1/d1-mat2/d2)^2)

 Correct me if I'm wrong though - aren't matrices of 48x times 48k
 larger than what R can handle at present?

 HTH

 Peter

 hmmm I didn't know that the limitation of R was below this value. I
 found this error message:

 Error in matrix(0, 48000, 48000) : too many elements specified

 but I thought it was a machine limitation (and I was asking for access
 to a better machine in my labs...). Thanks for clarifying it.

 Well, when Sarah gave me the answer for my problem, I got a new one :)
 Thank you, Sarah and Peter.

 So, is there any other way to trick R and allocate such large matrices?

 Best,

 Charles
 --
 Um axé! :)


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

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


[R] How to find matching columns in a matrix of lists?

2012-12-04 Thread Asis Hallab
Dear R users,

I have a matrix composed of lists:

m - matrix( list(), nrow=1, ncol=3 )
m[[ 1, 1 ]] - list(A, B)
m[[ 1, 2 ]] - list(A, C)
m[[ 1, 3 ]] - list(A, B)

and want to get the sub-matrix where cells contain B.
But

m[ , B %in% m[ 1, ], drop=F ]

as well as

m[ , B %in% m[ 1, ][], drop=F ]

return empty matrices.

Any ideas, hints and help will be very much appreciated!
Kind regards!
Josef

[[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] reformatting some data

2012-12-04 Thread Charles Determan Jr
Hello,

I am trying to reformat some data so that it is organized by group in the
columns.  The data currently looks like this:

   group X3.Hydroxybutyrate X3.Hydroxyisovalerate   ADP
347 4  4e-04 3e-04  5e-04
353 3  5e-04 3e-04  6e-04
359 4  4e-04 3e-04  6e-04
365 4  6e-04 3e-04  5e-04
371 4  5e-04 3e-04  7e-04
377 2  7e-04 4e-04  7e-04

I would like to reformat it so it is like this:

2  3   4
var1
var2
var3


I realize that there unequal numbers in each group but I would like to
none-the-less if possible.
Here is a subset of the data:

structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
c(4e-04,
5e-04, 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
365L, 371L, 377L), class = data.frame)

Any insight is truly appreciated,
Regards,
Charles

[[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] reformatting some data

2012-12-04 Thread Mark Lamias
It's not clear to me what it is you are attempting to do, as you switch from a 
very specific example to some general example with the vague terms var1 
var2, and var3.

It sounds like you might be trying to do something similar to what would be 
available in the shape package using the melt function.

Try ?melt.data.frame.

--Mark Lamias





 From: Charles Determan Jr deter...@umn.edu
To: r-help@r-project.org 
Sent: Tuesday, December 4, 2012 4:17 PM
Subject: [R] reformatting some data

Hello,

I am trying to reformat some data so that it is organized by group in the
columns.  The data currently looks like this:

       group X3.Hydroxybutyrate X3.Hydroxyisovalerate   ADP
347     4              4e-04                 3e-04                  5e-04
353     3              5e-04                 3e-04                  6e-04
359     4              4e-04                 3e-04                  6e-04
365     4              6e-04                 3e-04                  5e-04
371     4              5e-04                 3e-04                  7e-04
377     2              7e-04                 4e-04                  7e-04

I would like to reformat it so it is like this:

                2          3           4
var1
var2
var3


I realize that there unequal numbers in each group but I would like to
none-the-less if possible.
Here is a subset of the data:

structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
c(4e-04,
5e-04, 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
365L, 371L, 377L), class = data.frame)

Any insight is truly appreciated,
Regards,
Charles

    [[alternative HTML version deleted]]

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

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


Re: [R] reformatting some data

2012-12-04 Thread Mark Lamias
That should read the reshape package -- not the shape package.  

My apologies.






To: Charles Determan Jr deter...@umn.edu; r-help@r-project.org 
r-help@r-project.org 
Sent: Tuesday, December 4, 2012 4:36 PM
Subject: Re: [R] reformatting some data


It's not clear to me what it is you are attempting to do, as you switch from a 
very specific example to some general example with the vague terms var1 
var2, and var3.

It sounds like you might be trying to do something similar to what would be 
available in the shape package using the melt function.

Try ?melt.data.frame.

--Mark Lamias





 From: Charles Determan Jr deter...@umn.edu
To: r-help@r-project.org 
Sent: Tuesday, December 4, 2012 4:17 PM
Subject: [R] reformatting some data

Hello,

I am trying to reformat some data so that it is organized by group in the
columns.  The data currently looks like this:

       group X3.Hydroxybutyrate X3.Hydroxyisovalerate   ADP
347     4              4e-04                 3e-04                  5e-04
353     3              5e-04                 3e-04                  6e-04
359     4              4e-04                 3e-04                  6e-04
365     4              6e-04           
      3e-04                  5e-04
371     4              5e-04                 3e-04                  7e-04
377     2              7e-04                 4e-04                  7e-04

I would like to reformat it so it is like this:

                2          3           4
var1
var2
var3


I realize that there unequal numbers in each group but I would like to
none-the-less if possible.
Here is a subset of the data:

structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
c(4e-04,
5e-04,
 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
365L, 371L, 377L), class = data.frame)

Any insight is truly appreciated,
Regards,
Charles

    [[alternative HTML version deleted]]

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

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


Re: [R] monte carlo simulation on R

2012-12-04 Thread Suzen, Mehmet
There is a book on MC  with R:

Introducing Monte Carlo Methods with R by Robert/Casella:
http://www.springer.com/statistics/computational+statistics/book/978-1-4419-1575-7

On 4 December 2012 19:38, R. Michael Weylandt
michael.weyla...@gmail.com wrote:
 replicate(1000, sum(rnorm(50)^2-rchisq(50, 3)))

 Or you know, many other things...

 Michael

 On Tuesday, December 4, 2012, Adel ESSAFI wrote:

 Hello,

 How can I make a monte carlo simulation on R?

 Regards
 Adel


 --
 PhD candidate in Computer Science
 Address
 3 avenue lamine, cité ezzahra, Sousse 4000
 Tunisia
 tel: +216 97 246 706 (+33640302046 jusqu'au 15/6)
 fax: +216 71 391 166

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

2012-12-04 Thread arun


HI,

I just wonder whether your code worked or not.
set.seed(8)
mat1-matrix(sample(1:80,40,replace=TRUE),ncol=8)
set.seed(25)
mat2-matrix(sample(1:160,40,replace=TRUE),ncol=8)
#Since the dimensions are the same, 
m-1:5
n-1:8
sum1-0

for(i in 1:length(m)){
 for(j in 1:length(n)){
 sum1-sum1+(((mat1[i,j]/d1)-(mat2[i,j]/d2))^2)
 }
 }
sum1
#[1] 15192.89

#Sara's code:
sum((mat1/d1 - mat2/d2)^2)
#[1] 15192.89
A.K.




- Original Message -
From: Charles Novaes de Santana charles.sant...@gmail.com
To: r-help@r-project.org r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 2:27 PM
Subject: Re: [R] Large loops in R

Dear Michael,

Thank you for your answer.

I have 2 matrices. Each position of the matrices is a weight. And I
need to calculate the following sum of differences:

Considering:
mat1 and mat2 - two matrices (each of them 48000 x 48000).
d1 and d2 - two constant values.

sum-0;
for(i in 1:nrows1){
                        for(j in 1:nrows2){
                                        sum-sum+ ( ( (mat1(i,j)/d1) -
(mat2(i,j)/d2) )^2 )
                                }
                        }
                }

I was wondering if there is a better way to do this sum.

Thank you for your attention!

Best,

Charles

On Tue, Dec 4, 2012 at 7:54 PM, R. Michael Weylandt
michael.weyla...@gmail.com michael.weyla...@gmail.com wrote:


 On Dec 4, 2012, at 6:47 PM, Charles Novaes de Santana 
 charles.sant...@gmail.com wrote:

 Dear all,

 I need to access data from a large matrix (48000 x 48000) and to do it
 I am trying to run two loops using for command. Surely it is been a
 very slow job.

 I heard that for is not the best option to perform large loops in R,
 but I don't really know what would be the best (fast) option. sapply?
 vapply? Could anyone help me with this issue, please?

 What exactly are you trying to do? It's likely doable with a few vectorized 
 operations.

 Michael



 Thank you very much for your attention and for any help!

 Best regards,

 Charles

 --
 Um axé! :)

 --
 Charles Novaes de Santana
 http://www.imedea.uib-csic.es/~charles
 PhD student - Global Change
 Laboratorio Internacional de Cambio Global
 Department of Global Change Research
 Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
 Calle Miquel Marques 21, 07190
 Esporles - Islas Baleares - España

 Office phone - +34 971 610 896
 Cell phone - +34 660 207 940

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



-- 
Um axé! :)

--
Charles Novaes de Santana
http://www.imedea.uib-csic.es/~charles
PhD student - Global Change
Laboratorio Internacional de Cambio Global
Department of Global Change Research
Instituto Mediterráneo de Estudios Avanzados(CSIC/UIB)
Calle Miquel Marques 21, 07190
Esporles - Islas Baleares - España

Office phone - +34 971 610 896
Cell phone - +34 660 207 940

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] computing marginal values based on multiple columns?

2012-12-04 Thread Simon
Thanks to both of you! I learned something new from both of your posts :-)

And you were both right, my example numbers were wrong! I accidentally
computed them based on 0.25 * mean instead of 0.2 * mean.

Thanks again!

Simon

On Wed, Dec 5, 2012 at 4:07 AM, arun smartpink...@yahoo.com wrote:

 HI,

 I am not sure the output you wanted is correct:

 
 sample1 sample2 sample3
 1  1.0 0 0.5
 

 because
 0.2*colMeans(x[,-4])
 sample1 sample2 sample3
 #  28.40   24.08   21.36


 This might help you:
 apply(x[-4],2,function(y) length(y[y 0.2*mean(y) 
 x$class==a])/length(x[x$class==a]))
 #sample1 sample2 sample3
   #  0.0 0.0 0.5
 A.K.



 - Original Message -
 From: Simon simonzm...@gmail.com
 To: r-help@r-project.org
 Cc:
 Sent: Tuesday, December 4, 2012 4:49 AM
 Subject: [R] computing marginal values based on multiple columns?

 Hello all,

 I have what feels like a simple problem, but I can't find an simple
 answer. Consider this data frame:

  x - data.frame(sample1=c(35,176,182,193,124),
 sample2=c(198,176,190,23,15), sample3=c(12,154,21,191,156),
 class=c('a','a','c','b','c'))

  x
   sample1 sample2 sample3 class
 1  35 198  12 a
 2 176 176 154 a
 3 182 190  21 c
 4 193  23 191 b
 5 124  15 156 c

 Now I wish to know: for each sample, for values  20% of the sample mean,
 what percentage of those are class a?

 I want to end up with a table like:

sample1 sample2 sample3
 1  1.0 0 0.5

 I can calculate this for an individual sample using this rather clumsy
 expression:

 length(which(x$sample1  mean(x$sample1)  x$class=='a')) /
 length(which(x$sample1  mean(x$sample1)))

 I'd normally propagate it across the data frame using apply, but I
 can't because it depends on more than one column.

 Any help much appreciated!

 Cheers,

 Simon

 [[alternative HTML version deleted]]

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



[[alternative HTML version deleted]]

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


Re: [R] list to matrix?

2012-12-04 Thread arun
Hi,
Try this:
list1-list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17), 
    c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25), 
    c(45, 19), c(5e+05, 16))


res-t(sapply(list1,function(x) x))
res
#    [,1] [,2]
 #[1,]  5  101
 #[2,] 10   46
 #[3,] 15   31
 #[4,] 20   17
 #[5,] 25   19
 #[6,] 30   11
 #[7,] 35   12
 #[8,] 40   25
 #[9,] 45   19
#[10,] 50   16

  is.matrix(res)
#[1] TRUE
#or
res1-sapply(list1,function(x) x)
 is.matrix(res1)
#[1] TRUE
A.K.


- Original Message -
From: Sam Steingold s...@gnu.org
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 3:09 PM
Subject: [R] list to matrix?

How do I convert a list to a matrix?

--8---cut here---start-8---
list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17), 
    c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25), 
    c(45, 19), c(5e+05, 16))
as.matrix(a)
      [,1]    
[1,] Numeric,2
[2,] Numeric,2
[3,] Numeric,2
[4,] Numeric,2
[5,] Numeric,2
[6,] Numeric,2
[7,] Numeric,2
[8,] Numeric,2
[9,] Numeric,2
--8---cut here---end---8---

thanks!

-- 
Sam Steingold (http://sds.podval.org/) on Ubuntu 12.04 (precise) X 11.0.11103000
http://www.childpsy.net/ http://palestinefacts.org http://dhimmi.com
http://jihadwatch.org http://www.PetitionOnline.com/tap12009/ http://memri.org
Rhinoceros has poor vision, but, due to his size, it's not his problem.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] How do I get internal nodes of dendograms produced by R?

2012-12-04 Thread sagnik ray choudhury
I am using R for hierarchical clustering of a number of documents.

I have a distance matrix on which I have applied hclust method. When I plot
the result of hclust method, I can see the dendogram plotted. What I need
now is the dendogram stored as a tree in a data structure. My goal is to
automatically label all internal nodes. For that, I need to know, which
leaf nodes make a first level cluster, and which first level clusters make
a second level cluster and so on. Is there a way in R to get this
information?



-- 
Sagnik Ray Choudhury
sag...@psu.edu

[[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] How to find matching columns in a matrix of lists?

2012-12-04 Thread Asis Hallab
Hello,

m[ , sapply(1:ncol(m), function(j) sapply(B, `%in%`, m[[1 , j]])), drop=F
 ]


It indeed does.

Thank you very much!

[[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] reformatting some data

2012-12-04 Thread arun
Hi,

Not sure whether this helps:
library(reshape2)

dat1-structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
 c(4e-04,
 5e-04, 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
 3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
 5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
 X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
 365L, 371L, 377L), class = data.frame)
datM-melt(dat1,id.var=group)
dcast(datM,variable~group,length)
#or
dcast(datM,variable~group,mean)
#   variable 2 3    4
#1    X3.Hydroxybutyrate 7e-04 5e-04 0.000475
#2 X3.Hydroxyisovalerate 4e-04 3e-04 0.000300
#3   ADP 7e-04 6e-04 0.000575
A.K




- Original Message -
From: Charles Determan Jr deter...@umn.edu
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 4:17 PM
Subject: [R] reformatting some data

Hello,

I am trying to reformat some data so that it is organized by group in the
columns.  The data currently looks like this:

       group X3.Hydroxybutyrate X3.Hydroxyisovalerate   ADP
347     4              4e-04                 3e-04                  5e-04
353     3              5e-04                 3e-04                  6e-04
359     4              4e-04                 3e-04                  6e-04
365     4              6e-04                 3e-04                  5e-04
371     4              5e-04                 3e-04                  7e-04
377     2              7e-04                 4e-04                  7e-04

I would like to reformat it so it is like this:

                2          3           4
var1
var2
var3


I realize that there unequal numbers in each group but I would like to
none-the-less if possible.
Here is a subset of the data:

structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
c(4e-04,
5e-04, 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
365L, 371L, 377L), class = data.frame)

Any insight is truly appreciated,
Regards,
Charles

    [[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] reformatting some data

2012-12-04 Thread arun
Hi,
You can also do this:
dat1-structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
 c(4e-04,
 5e-04, 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
 3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
 5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
 X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
 365L, 371L, 377L), class = data.frame)

datM-melt(dat1,id.var=group)

xtabs(value~variable+group,data=datM)
#   group
#variable 2  3  4
 # X3.Hydroxybutyrate    0.0007 0.0005 0.0019
 # X3.Hydroxyisovalerate 0.0004 0.0003 0.0012
 # ADP   0.0007 0.0006 0.0023
A.K.



- Original Message -
From: Charles Determan Jr deter...@umn.edu
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 4:17 PM
Subject: [R] reformatting some data

Hello,

I am trying to reformat some data so that it is organized by group in the
columns.  The data currently looks like this:

       group X3.Hydroxybutyrate X3.Hydroxyisovalerate   ADP
347     4              4e-04                 3e-04                  5e-04
353     3              5e-04                 3e-04                  6e-04
359     4              4e-04                 3e-04                  6e-04
365     4              6e-04                 3e-04                  5e-04
371     4              5e-04                 3e-04                  7e-04
377     2              7e-04                 4e-04                  7e-04

I would like to reformat it so it is like this:

                2          3           4
var1
var2
var3


I realize that there unequal numbers in each group but I would like to
none-the-less if possible.
Here is a subset of the data:

structure(list(group = c(4L, 3L, 4L, 4L, 4L, 2L), X3.Hydroxybutyrate =
c(4e-04,
5e-04, 4e-04, 6e-04, 5e-04, 7e-04), X3.Hydroxyisovalerate = c(3e-04,
3e-04, 3e-04, 3e-04, 3e-04, 4e-04), ADP = c(5e-04, 6e-04, 6e-04,
5e-04, 7e-04, 7e-04)), .Names = c(group, X3.Hydroxybutyrate,
X3.Hydroxyisovalerate, ADP), row.names = c(347L, 353L, 359L,
365L, 371L, 377L), class = data.frame)

Any insight is truly appreciated,
Regards,
Charles

    [[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] How to find matching columns in a matrix of lists?

2012-12-04 Thread arun
HI,

Does this work for you?
 mapply(function(x) x==B,m)
  [,1]  [,2]  [,3]
#[1,] FALSE FALSE FALSE
#[2,]  TRUE FALSE  TRUE
A.K.




- Original Message -
From: Asis Hallab asis.hal...@gmail.com
To: r-help@r-project.org
Cc: 
Sent: Tuesday, December 4, 2012 4:16 PM
Subject: [R] How to find matching columns in a matrix of lists?

Dear R users,

I have a matrix composed of lists:

m - matrix( list(), nrow=1, ncol=3 )
m[[ 1, 1 ]] - list(A, B)
m[[ 1, 2 ]] - list(A, C)
m[[ 1, 3 ]] - list(A, B)

and want to get the sub-matrix where cells contain B.
But

m[ , B %in% m[ 1, ], drop=F ]

as well as

m[ , B %in% m[ 1, ][], drop=F ]

return empty matrices.

Any ideas, hints and help will be very much appreciated!
Kind regards!
Josef

    [[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] How to find matching columns in a matrix of lists?

2012-12-04 Thread Rui Barradas

Hello,

Try the following.

m[ , sapply(1:ncol(m), function(j) sapply(B, `%in%`, m[[1 , j]])), 
drop=F ]



Hope this helps,

Rui Barradas
Em 04-12-2012 21:16, Asis Hallab escreveu:

Dear R users,

I have a matrix composed of lists:

m - matrix( list(), nrow=1, ncol=3 )
m[[ 1, 1 ]] - list(A, B)
m[[ 1, 2 ]] - list(A, C)
m[[ 1, 3 ]] - list(A, B)

and want to get the sub-matrix where cells contain B.
But

m[ , B %in% m[ 1, ], drop=F ]

as well as

m[ , B %in% m[ 1, ][], drop=F ]

return empty matrices.

Any ideas, hints and help will be very much appreciated!
Kind regards!
Josef

[[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] list to matrix?

2012-12-04 Thread R. Michael Weylandt
On Tue, Dec 4, 2012 at 8:17 PM, arun smartpink...@yahoo.com wrote:
 Hi,
 Try this:
 list1-list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17),
 c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25),
 c(45, 19), c(5e+05, 16))


 res-t(sapply(list1,function(x) x))

Bah Humbug! (In Christmas cheer)

No need for all this (see solutions including mine already given) --
but even without those, this is silly. An identity map is a real waste
if you just want the simplification bit of sapply() -- you'd be much
better just using simplify2array()

Michael

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] latticeExtra tileplot question - tiles are not all the same size, need help.

2012-12-04 Thread Ryan Flaherty
Hello,

I have been creating many tileplots to try and illustrate the relative
abundance of fish through space and time. My issue is that the tiles that
border the plot are smaller than those in the center of the plot. In the
example I've provided the effect is pretty minor (I'm hoping this will be an
adequate example as I had the code already created/data uploaded). However,
I have other plots that have fewer space-time strata (fewer tiles) and the
effect is very noticeable. 

Code:
##You will need to load the following packages: latticeExtra, gridExtra and
RCurl
#Upload data


cwtb_csv-getURL(https://docs.google.com/spreadsheet/pub?key=0AjzYZNH9Dw9qd
Fc5c1FXZ19uYXE3U1QwU1MxVkR4dGcsingle=truegid=0output=csv,ssl.verifypeer
= FALSE)
cwtb-read.csv(textConnection(cwtb_csv),header=T)
gsi_csv-getURL(https://docs.google.com/spreadsheet/pub?key=0AjzYZNH9Dw9qdD
NTN2djUEl3UVdhZ0t4ZXZrdHpPc2coutput=csv,ssl.verifypeer=FALSE)
gsi_s-read.csv(textConnection(gsi_csv),header=T)

##Custom color ramp...I didn't like the ones that are provided in the
package, this is similar to jet colors in matlab
jet.colors -
  colorRampPalette(c(#7F, blue, #007FFF, cyan,
 #7FFF7F, yellow, #FF7F00, red, #7F))

#Tileplots comparing how CWT recoveries from brood year 1974 -1977 and GSI
data 
#from 2010 - 2011 represent the relative abundance of Rogue river chinook
salmon through space and time

R_g_s-tileplot(Rogue~month_num*-area_num,gsi_s,col.regions=jet.colors(256),
scales=list(x=list(at=5:9, labels=c('May'
,'June','July','August','September')), 
 y=list(at=(-1):-8,labels=c('T','N','C','B','KC','FB','SF','MO'))),
main = list(label=Relative abundance of Rogue River stock (WCGSI
2010-2011),cex=0.75) ,
xlab= Month,
ylab= Management Area,
border = black,
panel = function(...){
panel.fill(black)
panel.voronoi(...)
})

R_c_b-tileplot(Rogue~month_num*-area_num,cwtb,col.regions=jet.colors(256),
scales=list(x=list(at=5:9, labels=c('May'
,'June','July','August','September')), 
  y=list(at=(-1):-8,labels=c('T','N','C','B','KC','FB','SF','MO'))),
main = list(label=Relative abundance of Rogue River stock (CWT 1977
- 1983),cex=0.75) ,
xlab= Month,
ylab= Management Area,
border = black,
panel = function(...){
panel.fill(black)
panel.voronoi(...)
})

#plot the two tileplots side by side

grid.arrange(R_c_b,R_g_s,ncol=2)

Thank you in advance for any advice regarding how to remedy this issue.

Regards,


Ryan

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

2012-12-04 Thread William Dunlap
 No need for all this (see solutions including mine already given) --
 but even without those, this is silly. An identity map is a real waste
 if you just want the simplification bit of sapply() -- you'd be much
 better just using simplify2array()

You are right that simplify2array(p) does everything that sapply(p,function(x)x)
does and does it more quickly and that matrix(unlist(p),nrow=2,byrow=TRUE)
is much faster than either:

   p - lapply(1:1e6, function(i)c(i, log2(i)))
   system.time(z1 - t(sapply(p,function(x)x)))
 user  system elapsed 
  1.2 0.0 1.2 
   system.time(z2 - t(simplify2array(p)))
 user  system elapsed 
  0.910.000.90 
   system.time(z3 - matrix(unlist(p), ncol=2, byrow=TRUE))
 user  system elapsed 
 0.040.000.04

You can also use vapply instead of sapply - it requires that you supply
the expect shape and type of FUN's output so it is doesn't have to figure
this out from looking at all the outputs of FUN: 
   system.time(z4 - t(vapply(p,FUN=function(x)x,FUN.VALUE=numeric(2
 user  system elapsed 
 0.560.000.56 

An advantage of vapply is that it stops in its tracks if FUN returns
an unexpected value.  sapply() and simplify2array() will silently give
you an unexpected result (a single column matrix of mode list 
instead of a vector of numbers)  and matrix(unlist()...) gives you a
warning if you are lucky.
   pBad - p ; pBad[[length(pBad)/2]] - 666
   system.time(zBad1 - t(sapply(pBad,function(x)x)))
 user  system elapsed 
 1.750.001.75 
   zBad1[,49:51] # not what we wanted
  [[1]]
  [1] 49.0 18.93157

  [[2]]
  [1] 666
  
  [[3]]
  [1] 51.0 18.93157
   system.time(zBad2 - t(simplify2array(pBad)))
 user  system elapsed 
  0.5 0.0 0.5 
   system.time(zBad3 - matrix(unlist(pBad), ncol=2, byrow=TRUE))
 user  system elapsed 
 0.030.000.03 
  Warning message:
  In matrix(unlist(pBad), ncol = 2, byrow = TRUE) :
data length [199] is not a sub-multiple or multiple of the number of 
rows [100]
   # no warning if length(unlist(p)) were even
   system.time(zBad4 - t(vapply(pBad,function(x)x,numeric(2
  Error in vapply(pBad, function(x) x, numeric(2)) : 
values must be length 2,
   but FUN(X[[50]]) result is length 1
  Timing stopped at: 0.29 0 0.28

Which of the latter two methods you choose depends on how likely errors
in the data are.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
 Behalf
 Of R. Michael Weylandt
 Sent: Tuesday, December 04, 2012 2:59 PM
 To: arun
 Cc: R help; s...@gnu.org
 Subject: Re: [R] list to matrix?
 
 On Tue, Dec 4, 2012 at 8:17 PM, arun smartpink...@yahoo.com wrote:
  Hi,
  Try this:
  list1-list(c(5, 101), c(1e+05, 46), c(15, 31), c(2e+05, 17),
  c(25, 19), c(3e+05, 11), c(35, 12), c(4e+05, 25),
  c(45, 19), c(5e+05, 16))
 
 
  res-t(sapply(list1,function(x) x))
 
 Bah Humbug! (In Christmas cheer)
 
 No need for all this (see solutions including mine already given) --
 but even without those, this is silly. An identity map is a real waste
 if you just want the simplification bit of sapply() -- you'd be much
 better just using simplify2array()
 
 Michael
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] partial analisys of a time series

2012-12-04 Thread Antonio Silva
Hi all

Thanks for the attention and answers. I learned a lot I now I can go on my
work. I also tryed to you de command window().

I thought it would be possible to select one column of an ts object, like
we can do with a data.frame (plot(data[,2],data[,3]), to work. But as I saw
we need to extract the values and create another ts.

Thanks very much. All the best.

Antonio



2012/12/4 arun smartpink...@yahoo.com

 Hi,
 If the frequency is 1, the error message will  be gone.
 For e.g.

 birthstimeseriesJanFeb-subset(birthstimeseries,cycle(birthstimeseries)==c(1,2))

  
 birthstimeseriesJanFeb1-ts(birthstimeseriesJanFeb,frequency=2,start=c(1946,1))
  plot.ts(birthstimeseriesJanFeb1)
  birthstimeseriesJanFebHW-HoltWinters(birthstimeseriesJanFeb1)
  plot(birthstimeseriesJanFebHW)
 A.K.




 - Original Message -
 From: Antonio Silva aolinto@gmail.com
 To: PIKAL Petr petr.pi...@precheza.cz
 Cc: R-help@r-project.org R-help@r-project.org
 Sent: Tuesday, December 4, 2012 5:50 AM
 Subject: Re: [R] partial analisys of a time series

 Thanks Petr

 I thought there might be an equivalent for birthstimeseries[,1] if it were
 a dataframe, but split function sounds great.

 I could not reproduce the second line of your suggestion l.blist -
 lapply(blist, HoltWinters). I receive the message: Error in
 decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) : time
 series has no or less than 2 periods

 What could be going wrong?

 Best regards

 Antonio


 2012/12/4 PIKAL Petr petr.pi...@precheza.cz

  Hi
 
   -Original Message-
   From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
   project.org] On Behalf Of Antonio Silva
   Sent: Tuesday, December 04, 2012 10:26 AM
   To: R-help@r-project.org
   Subject: [R] partial analisys of a time series
  
   Dear list members
  
   I want to analyze separately the months of a time series. In other
   words, I want to plot and fit models for each month separately.
  
   Taking the example of
   http://a-little-book-of-r-for-time-
   series.readthedocs.org/en/latest/src/timeseries.html
  
   births - scan(http://robjhyndman.com/tsdldata/data/nybirths.dat;)
   birthstimeseries - ts(births, frequency=12, start=c(1946,1))
   birthstimeseries
   plot.ts(birthstimeseries)
   birthstimeseriesHW - HoltWinters(birthstimeseries)
   plot(birthstimeseriesHW)
  
   How to proceed the plotting and HoltWinters smoothing considereing only
   Januarys, Februarys, etc. separately.
 
  Split your data by months to a list, use lapply.
 
  using zoo package
 
  blist -split(birthstimeseries, months(as.Date(birthstimeseries)))
  l.blist - lapply(blist, HoltWinters)
 
  Regards
  Petr
 
 
 
 
 
  
   Thanks in advance.
  
   Antonio Olinto
  
 [[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.
 



 --
 Antônio Olinto Ávila da Silva
 Biólogo / Oceanógrafo
 Instituto de Pesca (Fisheries Institute)
 São Paulo, Brasil

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




-- 
Antônio Olinto Ávila da Silva
Biólogo / Oceanógrafo
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

[[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] control point size of superscript when labeling axes with title()

2012-12-04 Thread David Winsemius

On Dec 4, 2012, at 11:05 AM, Chris Solomon wrote:

 Hi-
 
 A journal has asked me to make all of my text annotations on a figure at 
 10-point size. For the most part this is easy, e.g. by creating figures with:
   pdf(..., family='Times', pointsize=10)
 
 But where I have superscripts (or subscripts) in axis labels, the default 
 seems to be to shrink the superscripted text slightly. For example this code:
   title(ylab=expression(paste('Respiration (mg  ',O[2],'  ',L^-1,' 
 ',d^-1,')',sep=' ')),outer=T,line=0.3)
 produces superscripted numbers at approximately 7 point.
 
 I have been poking around for a solution but not having much luck. The 
 textstyle() function keeps superscripted text at original size if you use it 
 within a text() call, but I can't figure out an equivalent solution within a 
 title() call.

It is not the fact that you are using title, but rather (I think) the fact that 
numbers are handled differently than text in plotmath. You are also 
unnecessarily using 'paste' and incorrectly including 'sep' argument in 
plotmath-paste, which does not honor that argument. 

Regular R-paste != plotmath-paste

Quoting the exponents allows textstyle to do its job properly. 

 plot(1,1, ylab=)
 title(ylab = 
 expression(Respiration~(*mg~~O[2]~~L^textstyle(-1)~d^textstyle(-1)*')') 
 ,line= 1.5)

   [[alternative HTML version deleted]]
 
You might want to read what the Posting Guide says about HTML posting.

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

David Winsemius, MD
Alameda, CA, 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] How to find matching columns in a matrix of lists?

2012-12-04 Thread Asis Hallab
Hi R users following this thread!

I evaluated both solutions given by Rui and Arun.
Both work very well.
Arun's is a little faster.

I did the following time measurements on a large matrix with 2 rows,
labelled InterPro and GO, and 88664 columns, labelled with protein IDs.
I was interested in selecting those columns (proteins) that share a certain
GO annotation for molecular function.

So here is the time evaluation:
* Rui's method:
system.time(annos[ , sapply( annos[ GO, ], function(x) { GO:0005634
%in% x } ), drop=F ] )
   user  system elapsed
  1.068   0.012   1.079
or
system.time(annos[ , unlist( lapply( annos[ GO, ], function(x) {
GO:0005634 %in% x } ) ), drop=F ])
   user  system elapsed
  0.876   0.000   0.880

* Arun's method:
system.time( annos[ , mapply( function(x){any(x==GO:0005634)}, annos[
GO, ] ), drop=F ] )
   user  system elapsed
  0.808   0.012   0.826

Thank you Rui and Arun very much for your help and everyone else for your
attention!
Kind regards!

2012/12/5 arun smartpink...@yahoo.com

 HI,
 Just tweaking my code also gives the same result:
 m[,mapply(function(x) any(x==B),m),drop=F]
 # [,1]   [,2]
 #[1,] List,2 List,2
 A.K.



 - Original Message -
 From: Asis Hallab asis.hal...@gmail.com
 To: Rui Barradas ruipbarra...@sapo.pt; r-help@r-project.org
 Cc:
 Sent: Tuesday, December 4, 2012 5:26 PM
 Subject: Re: [R] How to find matching columns in a matrix of lists?

 Hello,

 m[ , sapply(1:ncol(m), function(j) sapply(B, `%in%`, m[[1 , j]])), drop=F
  ]
 

 It indeed does.

 Thank you very much!

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




-- 
Asis Hallab
Rothehausstr. 6 - 12
50823 Köln

Skype: asis.hallab.cgn
Fest (Köln) 42346046
Mobil  (O2) 0176 63370211
Fax 01212 - 5 - 30697106

[[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] control point size of superscript when labeling axes with title()

2012-12-04 Thread Chris Solomon
Fantastic - that does the trick, David, thanks for setting me straight!

Chris

-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net] 
Sent: Tuesday, December 04, 2012 8:19 PM
To: Chris Solomon
Cc: r-help@r-project.org
Subject: Re: [R] control point size of superscript when labeling axes with 
title()


On Dec 4, 2012, at 11:05 AM, Chris Solomon wrote:

 Hi-
 
 A journal has asked me to make all of my text annotations on a figure at 
 10-point size. For the most part this is easy, e.g. by creating figures with:
   pdf(..., family='Times', pointsize=10)
 
 But where I have superscripts (or subscripts) in axis labels, the default 
 seems to be to shrink the superscripted text slightly. For example this code:
   title(ylab=expression(paste('Respiration (mg  ',O[2],'  ',L^-1,' 
 ',d^-1,')',sep=' ')),outer=T,line=0.3) produces superscripted numbers at 
 approximately 7 point.
 
 I have been poking around for a solution but not having much luck. The 
 textstyle() function keeps superscripted text at original size if you use it 
 within a text() call, but I can't figure out an equivalent solution within a 
 title() call.

It is not the fact that you are using title, but rather (I think) the fact that 
numbers are handled differently than text in plotmath. You are also 
unnecessarily using 'paste' and incorrectly including 'sep' argument in 
plotmath-paste, which does not honor that argument. 

Regular R-paste != plotmath-paste

Quoting the exponents allows textstyle to do its job properly. 

 plot(1,1, ylab=)
 title(ylab = 
 expression(Respiration~(*mg~~O[2]~~L^textstyle(-1)~d^textstyle(-1
 )*')') ,line= 1.5)

   [[alternative HTML version deleted]]
 
You might want to read what the Posting Guide says about HTML posting.

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

David Winsemius, MD
Alameda, CA, 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.


  1   2   >