Re: [R] How to remove rows based on frequency of factor and then difference date scores

2010-08-25 Thread Chris Beeley
Many thanks to you both. I have now filed away for future reference the 2 
factor tapply as well as the extremely useful looking plyr library. And the 
code worked beautifully :-)



On 24 Aug 2010, at 19:47, Abhijit Dasgupta, PhD aikidasgu...@gmail.com 
wrote:

 The paste-y argument is my usual trick in these situations. I forget that 
 tapply can take multiple ordering arguments :)
 
 Abhijit
 
 On 08/24/2010 02:17 PM, David Winsemius wrote:
 
 On Aug 24, 2010, at 1:59 PM, Abhijit Dasgupta, PhD wrote:
 
 The only problem with this is that Chris's unique individuals are a 
 combination of Type and ID, as I understand it. So Type=A, ID=1 is a 
 different individual from Type=B,ID=1. So we need to create a unique 
 identifier per person, simplistically by uniqueID=paste(Type, ID, sep=''). 
 Then, using this new identifier, everything follows.
 
 I see your point. I agree that a tapply method should present both factors 
 in the indices argument.
 
  new.df - txt.df[ -which( txt.df$nn =1), ]
  new.df - new.df[ with(new.df, order(Type, ID) ), ]  # and possibly needs 
  to be ordered?
  new.df$diffdays - unlist( tapply(new.df$dt2, list(new.df$ID, 
  new.df$Type), function(x) x[1] -x) )
  new.df
  Type ID   Date Valuedt2 nn diffdays
 1A  1 16/09/2020 8 2020-09-16  30
 2A  1 23/09/2010 9 2010-09-23  3 3646
 4B  1  13/5/2010 6 2010-05-13  30
 
 But do not agree that you need, in this case at least, to create a paste()-y 
 index. Agreed, however, such a construction can be useful in other 
 situations.
 
 
 
 -- 
 
 Abhijit Dasgupta, PhD
 Director and Principal Statistician
 ARAASTAT
 Ph: 301.385.3067
 E: adasgu...@araastat.com
 W: http://www.araastat.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] how do R calculates the number of intervals between tick-marks

2010-08-25 Thread Antonio Olinto

Hello,

I want to know how do R calculates the number of intervals between  
tick-marks in the y axis in a plot.


I'm making a three y-axes plot and this information would help me a lot.

Thanks in advance.

Antonio Olinto




Webmail - iBCMG Internet
http://www.ibcmg.com.br

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 remove all objects except a few specified objects?

2010-08-25 Thread Cheng Peng

Thanks Josh and Karl. function dnrm() works well for my purpose. 
-- 
View this message in context: 
http://r.789695.n4.nabble.com/How-to-remove-all-objects-except-a-few-specified-objects-tp2335651p2337398.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] Documenting S4 Methods

2010-08-25 Thread Dario Strbenac
I'm in the process of converting some S3 methods to S4 methods.

I have this function :

setGeneric(enrichmentCalc, function(rs, organism, seqLen, 
...){standardGeneric(enrichmentCalc)})

setMethod(enrichmentCalc, c(GenomeDataList, BSgenome), function(rs, 
organism, seqLen, ...) {
 ...  ... ...
})

setMethod(enrichmentCalc, c(GenomeData, BSgenome), function(rs, organism, 
seqLen=NULL, do.warn=FALSE) {
   ...... ...
})

setMethod(enrichmentCalc, c(GRanges, BSgenome), function(rs, organism, 
seqLen=NULL) {
   ...... ...
}

and a part of my Rd file is :

\name{enrichmentCalc}
\docType{methods}
\alias{enrichmentCalc,GenomeDataList,BSgenome-method}
\alias{enrichmentCalc,GenomeData,BSgenome-method}
\alias{enrichmentCalc,GRanges,BSgenome-method}
......   ...
\usage{
  enrichmentCalc(rs, organism, seqLen, ...)
  enrichmentCalc(rs, organism, seqLen=NULL, do.warn=FALSE)
  enrichmentCalc(rs, organism, seqLen=NULL)
}
.........

Can anyone suggest why I'm seeing this error :

* checking for code/documentation mismatches ... WARNING
Codoc mismatches from documentation object 'enrichmentCalc':
enrichmentCalc
  Code: function(rs, organism, seqLen, ...)
  Docs: function(rs, organism, seqLen = NULL, do.warn = FALSE)
  Argument names in code not in docs:
...
  Argument names in docs not in code:
do.warn
  Mismatches in argument names:
Position: 4 Code: ... Docs: do.warn
  Mismatches in argument default values:
Name: 'seqLen' Code:  Docs: NULL
enrichmentCalc
  Code: function(rs, organism, seqLen, ...)
  Docs: function(rs, organism, seqLen = NULL)
  Argument names in code not in docs:
...
  Mismatches in argument default values:
Name: 'seqLen' Code:  Docs: NULL

* checking Rd \usage sections ... WARNING
Objects in \usage without \alias in documentation object 'enrichmentCalc':
  enrichmentCalc

Also, what is the difference between

...  ...  ...
\docType{methods}
...  ...  ...
\alias{methodName,class-method}
...  ...  ...
\usage{methodName(arg1)}
...  ...  ...

and

...  ...  ...
\alias{methodName,class-method}
...  ...  ...
\usage
{
\S4method{methodName}{class}(arg1)
}
...  ...  ...

I've seen both ways used for S4 methods and don't know what is the underlying 
difference.

I haven't been able to find any good tutorials for the new S4 architecture 
(written post 2006), so I'm not sure where to start with S4.

Thanks,
   Dario.

--
Dario Strbenac
Research Assistant
Cancer Epigenetics
Garvan Institute of Medical Research
Darlinghurst NSW 2010
Australia

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

2010-08-25 Thread Jonathan Greenberg
Rhelpers:

I'm trying to make a barchart of a 2-group dataset
(barchart(x~y,data=data,groups=z,horizontal=FALSE)).  My problem is
that I can't, for the life of me, seem to get rid of the inter-bar
space -- box.ratio set to 1 doesn't do much.  Any ideas?  I'd
ideally want zero space between the bars.  Thanks!

--j

-- 
Jonathan A. Greenberg, PhD
Assistant Project Scientist
Center for Spatial Technologies and Remote Sensing (CSTARS)
Department of Land, Air and Water Resources
University of California, Davis
One Shields Avenue
Davis, CA 95616
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307

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

2010-08-25 Thread david dav
Hi,
I would like to emphasize (zoom) the zone of a nomogram where the
probability are  0.01
(nomogram built with nomogram, Design).
As a consequence, I don't need to draw the part of the Total points
axis with score  60 equivalent in my case to a linear predictor  4.5

- As far as I know, this is not possible with the arguments of the function.
- Changing the code of the function is beyond my abilities
 -- can not even create a nomo.f function with the same body:
  body(nomo) - expression({
  conf.lp - match.arg(conf.lp) . rest of the function body
this does not even work
 -- I am not sure to find the piece of code responsible for
defining the axis

nomogram(logistic.fit, fun = plogis, fun.at = c(seq(0.01,1,by=.2)),
lmgp=.2, lp = T,maxscale = 100,   total.sep.page = T )

Thanks
David

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


Re: [R] Plot bar lines like excel

2010-08-25 Thread abotaha

Thanks a lot for the nice explanation. however, you mention that it will be
nicer if there is data. yes you are right it be very kind of you if you
include a simple example. where my data are:

x1-c(11.5,9.38,9.3,9.8,9,9.06,9.42,8.8,9.05,8.14,8.2,7.59,6.92,6.47,7.12,
7.47,6.81,8.41,9.64,9.62,9.2,8.92,8,7.6,7.6,7.19,6.94,6.91,7,7.25,6.6,7.18,7.78,
7.37,7.29,6.71,7.05,6.69,6.05,6.38,6.18,7.67,7.34,7.13,6.95,6.8,6.45,6.81,6.27,6.35)

x2-c(11.5,8.959535,9.728104,9.609262,8.755494,9.221545,8.244458,7.63944,7.92052,7.690449,7.853247,
7.239616,6.609007,5.92946,6.47822,5.906936,6.966004,8.630517,9.506582,9.57479,9.236638,9.581875,8.838992,8.368731,8.608884,7.998023,7.918123,7.832322,7.930177,7.479222,6.866978,7.454062,8.206185,7.344037,7.059774,6.547646,7.005803,6.623987,5.992691,6.313481,6.194613,6.224266,7.084932,6.568976,6.43866,5.70081,6.7593,6.6929,6.46806,5.964816)

regards, 



-- 
View this message in context: 
http://r.789695.n4.nabble.com/Plot-bar-lines-like-excel-tp2337089p2337394.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] Odp: Finding pairs

2010-08-25 Thread Petr PIKAL
Hi

without other details it is probably impossible to give you any reasonable 
advice. Do you have your data already in R? What is their form? Are they 
in 2 columns in data frame? How did you get them paired?

So without some more information probably nobody will invest his time as 
it seems no trivial to me.

Regards
Petr

r-help-boun...@r-project.org napsal dne 24.08.2010 20:28:42:

 
 
 
 
 Dear R Helpers,
 
 
 I am a newbie and recently got introduced to R. I have a large database 
 containing the names of bank branch offices along-with other details. I 
am 
 into Operational Risk as envisaged by BASEL II Accord. 
 
 
 I am trying to express my problem and I am using only an indicative data 
which
 comes in coded format.
 
 
 
 
 A (branch)  B (controlled by)
 
 
 144   
 145  
 146   
 147   144 
 148   145 
 149   147
 151   146  
  ..  ...
  
 ..  ...
 
 
 where 144's etc are branch codes in a given city and B is subset of A.
 
 
 
 
 If a branch code appearing in A also appears in B (which is paired 
with 
 some otehr element of A e.g. 144 appearing in A, also appears in B and 
is 
 paired with 147 of A and likewise), then that means 144 is controlling 

 operations of bank office 147. Again, 147 itself appears again in B and 
is 
 paired with bank branch coded 149. Thus, 149 is controlled by 147 and 
147 is 
 controlled by 144. Likewise there are more than 700 hundred branch name 
codes available.
 
 
 My objective is to group them as follows -
 
 
 Bank Branch
 
 
 144  147149 
 
 
 145
 
 
 146   151  
 
 
 148
 .
 
 
 or even the following output will do.
 
 
 144
 147
 149
 
 
 145
 
 
 146
 151
 
 
 148
 151
 ..
 
 
 I understand I should be writing some R code to begin with which I had 
tried 
 also but as of now I am helpless. Please guide me.
 
 
 Mike
 
 
 
 
 
[[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 include trend (drift term) in arima.sim

2010-08-25 Thread StephenRichards

Thanks to Mark Leeds and Dennis Murphy for their suggestions.  The function
arima.sim() only simulates stationary series without a trend, so the best
approach appears to be to add the simulated stationary part to the trend as
follows:

  Temp - arima.sim(n=N.Years.Forecast, list(ar=AR.Coef, ma=MA.Coef,
intercept=Intercept), sd=SD)
  if (1 == D.) Simulated.Kappa - Trend + cumsum(Temp)
  if (2 == D.) Simulated.Kappa - Trend + cumsum(cumsum(Temp))

This appears to work well for d=1 in an ARIMA(p,d,q) model, but less well
for d=2, where the results appear to be unstable.  The problem is that if
the forecast starts with two or three biggish values of the differences then
these seem to take over the forecast of Kappa.
-- 
View this message in context: 
http://r.789695.n4.nabble.com/How-to-include-trend-drift-term-in-arima-sim-tp2331581p2337729.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] how to plot y-axis on the right of x-axis

2010-08-25 Thread Gavin Simpson
On Wed, 2010-08-25 at 07:12 +0800, elaine kuo wrote:
 Dear List,
 
 I have a richness data distributing across 20 N to 20 S latitude. (120 E-140
 E longitude).
 
 
 I would like to draw the richness in the north hemisphere and a regression
 line in the plot
 (x-axis: latitude, y-axis: richness in the north hemisphere).
 The above demand is done using plot.
 
 Then, south hemisphere richness and regression are required to be generated
 using
 the same y-axis above but an x-axis on the left side of the y-axis.
 (The higher latitude in the south hemisphere, the left it would move)

Elaine,

It is not very clear what it is that you want here. Two plot panels on
the same device like:

layout(1:2)
plot(1:100)
plot(1:100)
layout(1)

The second panel for southern species richness, do you mean you want the
plot to go like this:

plot(1:100, xlim = c(100,1))

i.e. back to front to the normal axis limits? If so, you'll need to
produce a regression model for southern richness just as you would for
northern richness, but then predict from it at a sequence of latitude
values that covers the range of your southern hemisphere observations,
then plot these predictions versus the sequence of latitudes using
lines().

For example, say you have:

set.seed(123)
rich - data.frame(richness = rpois(100, lambda = 5),
   latitude = sample(seq(0, 90, length = 1000),
  100))
plot(richness ~ latitude, data = rich, xlim = c(90,0))

## ignoring that this is a count and might be more meaningfully
## modelled using a Poisson GLM
mod - lm(richness ~ latitude, data = rich)

## new prediction data, 100 equally spaced latitudes
pdat - data.frame(latitude = seq(0, 90, length = 100))
pdat - within(pdat, richness - predict(mod, newdata = pdat))
lines(richness ~ latitude, data = pdat, col = red)

If that doesn't help, we'll need more info.

HTH

G

 Please kindly share how to design the south plot and regression line for
 richness.
 Also, please advise if any more info is in need.
 
 Elaine
 
   [[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.

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

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


[R] multi day R run crashed - is there a log

2010-08-25 Thread Martin Tomko

Dear all,
I am using an R 2.10 installation on a Windows 203 server that I have no 
control over. After a multi-day run I found that it was 
terminated/crashed. Is there any log kept by R where I could see whether 
something/what happened? The same process has been run beofre on a 
smaller dataset (also at least a day of computing) without problems.


Thanks
Martin

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


Re: [R] Draw a perpendicular line?

2010-08-25 Thread Benno Pütz
Maybe

perp.slope = -1/slope
abline(cy - cx*perp.slope, perp.slope)

where cx, cy are x- and y-coordinate of C, resp., and slope the slope you 
calculated for the line through A and B


Am 24.Aug.2010 um 0:04 schrieb CZ:

 
 Hi,
 
 I am trying to draw a perpendicular line from a point to two points.
 Mathematically I know how to do it, but to program it, I encounter some
 problem and hope can get help.  Thanks. 
 
 I have points, A, B and C.  I calculate the slope and intercept for line
 drawn between A and B.
 I am trying to check whether I can draw a perpendicular line from C to line
 AB and get the x,y value for the point D at the intersection. 
 
 Assume I get the slope of the perpendicular line, I will have my point (D)
 using variable x and y which is potentially on line AB.   My idea was using
 |AC|*|AC| = |AD|*|AD|+ |CD|*|CD|.  I don't know what function I may need to
 call to calculate the values for point D (uniroot?).
 
 Thank you. 
 
 
 
 -- 
 View this message in context: 
 http://r.789695.n4.nabble.com/Draw-a-perpendicular-line-tp2335882p2335882.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] Odp: Finding pairs

2010-08-25 Thread Mike Rhodes
Dear Mr Petr Pikal

I am extremely sorry for the manner I have raised the query. Actually that was 
my first post to this R forum and in fact even I was also bit confused while 
drafting the query, for which I really owe sorry to all for consuming the 
precious time. Perhaps I will try to redraft my query in a better way as 
follows. 

I have two datasets A and B containing the names of branch offices of a 
particular bank say XYZ plc bank. The XYZ bank has number of main branch 
offices (say Parent) and some small branch offices falling under the purview of 
these main branch offices (say Child).

The datalist A and B consists of these main branch office names as well as 
small branch office names. B is subset of A and these branch names are coded. 
Thus we have two datasets A and B as (again I am using only a
 portion of a large database just to have some idea)


A B
144  
145   
146   
147                  144    
148  145  
 
149  147
151  148



Now the branch 144 appears in A as well as in B and in B it is mapped with 147. 
This means branch 147 comes under the purview of main branch 144. Again 147 is 
controlling the branch 149 (since 147 also has appeared in B and is mapped with 
149 of A).

Similarly, branch 145 is controlling branch 148 which further controls 
operations of bank branch 151 and like wise.

So in the end I need an output something like -

Main Branch   Branch office1 Branch
 office2
144 147 149
145 148 151    
146 NA
  NA   
...
..

 
I understand again I am not able to put forward my query properly. But I must 
thank all of you for giving a patient reading to my query and for reverting 
back earlier. Thanks once again.

With warmest regards

Mike 


--- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:

From: Petr PIKAL petr.pi...@precheza.cz
Subject: Odp: [R] Finding
 pairs
To: Mike Rhodes mike_simpso...@yahoo.co.uk
Cc: r-help@r-project.org
Date: Wednesday, 25 August, 2010, 6:39

Hi

without other details it is probably impossible to give you any reasonable 
advice. Do you have your data already in R? What is their form? Are they 
in 2 columns in data frame? How did you get them paired?

So without some more information probably nobody will invest his time as 
it seems no trivial to me.

Regards
Petr

r-help-boun...@r-project.org napsal dne 24.08.2010 20:28:42:

 
 
 
 
 Dear R Helpers,
 
 
 I am a newbie and recently got introduced to R. I have a large database 
 containing the names of bank branch offices along-with other details. I 
am 
 into Operational
 Risk as envisaged by BASEL II Accord. 
 
 
 I am trying to express my problem and I am using only an indicative data 
which
 comes in coded format.
 
 
 
 
 A (branch)                      B (controlled by)
 
 
 144                   
 145                      
 146                   
 147                                       144 
 148                                       145 
 149       
                                147
 151                                       146  
  ..                                      ...
  
 ..                                      ...
 
 
 where 144's etc are branch codes in a given city and B is subset of A.
 
 
 
 
 If a branch code appearing in A also appears in B (which is paired 
with 
 some otehr element of A e.g. 144 appearing in A, also appears in B and 
is 
 paired with 147 of A and
 likewise), then that means 144 is controlling 

 operations of bank office 147. Again, 147 itself appears again in B and 
is 
 paired with bank branch coded 149. Thus, 149 is controlled by 147 and 
147 is 
 controlled by 144. Likewise there are more than 700 hundred branch name 
codes available.
 
 
 My objective is to group them as follows -
 
 
 Bank Branch
 
 
 144      147    149 
 
 
 145
 
 
 146       151  
 
 
 148
 .
 
 
 or even the following output will do.
 
 
 144
 147
 149
 
 
 145
 
 
 146
 151
 
 
 148
 151
 ..
 
 
 I understand I should be writing some R
 code to begin with which I had 
tried 
 also but as of now I am helpless. Please guide me.
 
 
 Mike
 
 
 
 
 
    [[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

2010-08-25 Thread Loos, Martin
Good afternoon!

It may seem trivial to some/most of You, but I found it difficult to properly 
include a C++-based .dll into a package that I want to build for usage in R. I 
read through the Writing R extensions...  R administration ... 
instructions, but it seems I did not grasp the bigger picture.

The way I figured out to use the .dll in my package finally worked for using 
that package from the R console, but it gives multiple error and warning 
messages when running the RCMD check pkg command in the shell.

What I did: I created the package skeleton via package.skeleton() including a 
function that makes both (1) a
library.dynam() call to a .dll and (2) a .C() call including the name of the 
function in that DLL. Next, I run RCMD INSTALL pkg to install the package and 
then created a /lib folder in the installed package containing the .dll file.
I tried alternatives, e.g., (a) containing a /src folder in may package 
skeleton with the .dll in it and/or (b) including zzz.R with a .First.lib 
function including the library call  but nothing worked out when running 
the package check.

help! thank You very much, MArtin

[[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] lmer() causes segfault

2010-08-25 Thread Bertolt Meyer

Ben Bolker bbolker at gmail.com writes:


Bertolt Meyer bmeyer at sozpsy.uzh.ch writes:



Hello lmer() - users,

A call to the lmer() function causes my installation of R (2.11.1 on
Mac OS X 10.5.8) to crash and I am trying to figure out the problem.


 [snip snip]


detach(package:nlme)
library(lme4)

mod1 - lmer(performance ~ time + (time | GroupID/StudentNumber),  
data

= dataset.long, na.action = na.omit)

However, this call results in a segfault:

*** caught segfault ***
address 0x154c3000, cause 'memory not mapped'

and a lengthy traceback. I can reproduce this error. It also occurs
when I don't load nlme before lme4. Can someone tell me what I am
doing wrong? Any help is greatly appreciated.


  This may well be a bug in lmer.  There have been a number of
fussy computational issues with the lme4 package on the Mac platform.


Ben, thanks for your reply. I tried to replicate this issue with a  
small clean data set on a windows machine. You can find the code for  
the data frame (100 observations from my data) at the end of this  
mail. Very simple: four test scores per student over time, and  
students are nested in groups. On my Windows installation, lmer()  
throws an error that does not seem to get caught on the Mac, resulting  
in the segfault:


library(lme4)
mlmoded1.lmer - lmer(Score ~ Time + (Time | GroupID/StudentID), data  
= test.data)


Error: length(f1) == length(f2) is not TRUE
Addditional Warnings:
1: In StudentID:GroupID :
  numeric expression has 100 elements: only first one is used
2: In StudentID:GroupID :
  numeric expression has 100 elements: only first one is used

It seems to me that I am committing a trivial error here and that I am  
too blind to see it. Any ideas?


Regards,
Bertolt


If it is at all possible, please (1) post the results of sessionInfo()
[which will in particular specify which version of lme4 you are  
using];
(2) possibly try this with the latest development version of lme4,  
from
R-forge, if that's feasible (it might be necessary to build the  
package

from source), and most importantly:

(3) create a reproducible (for others) example -- most easily by
posting your data on the web somewhere, but if that isn't possible
by simulating data similar to yours (if it doesn't happen with another
data set of similar structure, that's a clue -- it says it's some more
particular characteristic of your data that triggers the problem) and

(4) post to to *either* the R-sig-mac or the R-sig-mixed-models list,
where the post is more likely to come to the attention of those who
can help diagnose/fix ...

 good luck
   Ben Bolker




test.data - data.frame(c(17370, 17370, 17370, 17370, 17379, 17379,  
17379, 17379, 17387, 17387, 17387, 17387, 17391, 17391, 17391, 17391,  
17392, 17392, 17392, 17392, 17394, 17394, 17394, 17394, 17408, 17408,  
17408, 17408, 17419, 17419, 17419, 17419, 17429, 17429, 17429, 17429,  
17432, 17432, 17432, 17432, 17436, 17436, 17436, 17436, 17439, 17439,  
17439, 17439, 17470, 17470, 17470, 17470, 17220, 17220, 17220, 17220,  
17348, 17348, 17348, 17348, 17349, 17349, 17349, 17349, 17380, 17380,  
17380, 17380, 17398, 17398, 17398, 17398, 17400, 17400, 17400, 17400,  
17402, 17402, 17402, 17402, 17403, 17403, 17403, 17403, 17413, 17413,  
17413, 17413, 17416, 17416, 17416, 17416, 17420, 17420, 17420, 17420,  
17421, 17421, 17421, 17421), c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2,  
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,  
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), c(1, 2, 3, 4,  
1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3,  
4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2,  
3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1,  
2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4,  
1, 2, 3, 4), c(76.76, 81.83, 89.78, 92.82, 75.86, 81.84, 88.96, 92.28,  
75.28, 80.68, 88.62, 92.29, 76.60, 84.59, 92.03, 94.05, 75.57, 79.94,  
86.11, 90.25, 74.54, 81.42, 87.50, 90.71, 76.02, 83.68, 91.11, 94.14,  
76.31, 83.76, 90.44, 94.58, 72.29, 80.51, 86.09, 90.41, 74.99, 82.28,  
88.77, 92.26, 75.28, 81.92, 89.25, 92.64, 76.31, 83.93, 91.00, 94.60,  
76.31, 82.44, 90.57, 95.17, 76.94, 82.21, 83.81, 85.00, 79.96, 81.92,  
86.32, 90.05, 82.01, 84.81, 88.79, 93.10, 77.87, 82.94, 86.86, 90.31,  
77.87, 79.64, 85.66, 86.97, 79.35, 80.44, 84.26, 83.62, 79.06, 81.56,  
85.00, 87.43, 79.34, 81.47, 83.23, 86.86, 79.44, 80.37, 84.36, 89.11,  
78.77, 81.02, 81.60, 87.21, 75.75, 79.35, 80.38, 86.87, 76.04, 80.57,  
83.36, 86.31))


names(test.data) - c(StudentID, GroupID, Time, Score)

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

Re: [R] Draw a perpendicular line?

2010-08-25 Thread baptiste auguie
hi,

also, make sure you have set the aspect ratio to 1:1 when plotting (asp=1).

HTH,

baptiste

On 25 August 2010 10:20, Benno Pütz pu...@mpipsykl.mpg.de wrote:
 Maybe

 perp.slope = -1/slope
 abline(cy - cx*perp.slope, perp.slope)

 where cx, cy are x- and y-coordinate of C, resp., and slope the slope you 
 calculated for the line through A and B


 Am 24.Aug.2010 um 0:04 schrieb CZ:


 Hi,

 I am trying to draw a perpendicular line from a point to two points.
 Mathematically I know how to do it, but to program it, I encounter some
 problem and hope can get help.  Thanks.

 I have points, A, B and C.  I calculate the slope and intercept for line
 drawn between A and B.
 I am trying to check whether I can draw a perpendicular line from C to line
 AB and get the x,y value for the point D at the intersection.

 Assume I get the slope of the perpendicular line, I will have my point (D)
 using variable x and y which is potentially on line AB.   My idea was using
 |AC|*|AC| = |AD|*|AD|+ |CD|*|CD|.  I don't know what function I may need to
 call to calculate the values for point D (uniroot?).

 Thank you.



 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Draw-a-perpendicular-line-tp2335882p2335882.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.




-- 


Dr. Baptiste Auguié

Departamento de Química Física,
Universidade de Vigo,
Campus Universitario, 36310, Vigo, Spain

tel: +34 9868 18617
http://webs.uvigo.es/coloides

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

2010-08-25 Thread e-letter
On 24/08/2010, David Winsemius dwinsem...@comcast.net wrote:

 On Aug 24, 2010, at 9:37 AM, e-letter wrote:

 Readers,

 According to the documentation for the function 'plotmath' there is no
 apparent possibility to add the percent sign (%) to a plot function,

 Where did you see an assertion made???

Within R I entered the command:

?plotmath

Also accessed using:

help.start(browser=opera)

Navigating the web browser page:

packages
packages in /usr/lib/R/library
grdevices
plotmath

In the list headed 'syntax' and 'meaning' within the section 'details'.

 e.g.

 plot(a[,1]~b[,2],ylab=expression(x~%),xlab=expression(z))



 How to achieve this please?

 Read the plotmath helo page more carefully. The section immediatedly
 below the plotmath expressions points you to the use of the symbol()
 expression-function and to the points help page where generation of
 the available glyphs proceeds according to the advice on help(plotmath):

In my system the paragraph immediately after the list of features
(i.e. 'syntax','meaning') describes a note to TeX users. I cannot see
reference to 'symbol()'.

   TestChars - function(sign=1, font=1, ...)
 + {
 +if(font == 5) { sign - 1; r - c(32:126, 160:254)
 +} else if (l10n_info()$MBCS) r - 32:126 else r - 32:255
 +if (sign == -1) r - c(32:126, 160:255)
 +par(pty=s)
 +plot(c(-1,16), c(-1,16), type=n, xlab=, ylab=,
 + xaxs=i, yaxs=i)
 +grid(17, 17, lty=1)
 +for(i in r) try(points(i%%16, i%/%16, pch=sign*i, font=font,...))
 + }
   TestChars(font=5)

 Notice that the % sign is three characters to the right (i.e.
 higher) of the forall symbol that is produced by the example code

I can't see 'forall' in the code above.

 they offer. (The numbering proceeds from bottom up which confused me
 at first.)

What numbering?

The documentation makes reference to the command:

demo(plotmath)

I applied this command and could not see an instruction to produce the
percent (%) symbol.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Rotate x-axis label on log scale

2010-08-25 Thread Tim Elwell-Sutton
Yes, I do want the labels on the x-axis but the problem arises when the
y-axis is logarithmic presumably because the position labels on the x-axis
still needs to be defined in terms of x and y coordinates. 

The simplified examples below should make the problem clearer.
I can't see a way of altering the staxlab code to fix this but would be very
grateful for any suggestions.


#Example 1
#plotting on a linear scale the labels are rotated as desired
library(plotrix)
x - y - 1:4
labels= c('label1', 'label2', 'label3', 'label4')
plot(x,y, xaxt='n')
staxlab(side=1,at=1:4,labels=labels, srt=45)


#Example 2
#plotting on a y-axis log scale the labels don't appear
plot(x,y, xaxt='n', log='y')
staxlab(side=1,at=1:4,labels=labels, srt=45)


-Original Message-
From: David Winsemius [mailto:dwinsem...@comcast.net] 
Sent: Wednesday, August 25, 2010 11:46 AM
To: tesutton
Cc: 'Jim Lemon'; r-help@r-project.org
Subject: Re: [R] Rotate x-axis label on log scale


On Aug 24, 2010, at 11:37 PM, Tim Elwell-Sutton wrote:

 Hi Jim
 Thanks for this. The staxlab function seems very useful.  
 Unfortunately, the
 rotation option doesn't seem to work for me when the y-axis is on a  
 log
 scale.

What part of:

If srt is not NA, the labels will be rotated srt degrees and placed  
below the plot. This method will only place labels at the bottom.

... is unclear? You did say you wanted the rotation to be on the x- 
axis, did you not?

You could, of course, look at Lemon's code and hack it to do the y-axis:

else {
 xylim - par(usr)
 ypos - xylim[3] - ticklen * (xylim[4] - xylim[3])
 par(xpd = TRUE)
 text(at, ypos, labels, srt = srt, adj = 1, ...)
 par(xpd = FALSE)
 }


-- 
david.


 It will stagger the labels but not rotate them. There's no error
 message. On a linear axis the rotation works nicely. Any ideas?
 The example below works if you omit log='y' or srt=45

 Thanks very much
 Tim

 #Create plot with log-scale on the y-axis
 par(mar = c(7, 4, 4, 2) + 0.1)
 plot(1, type='n', bty='n',
xlab=,
ylab='Odds Ratio',
xlim= c(0.5,4.5),
ylim= c(0.75, 2),
cex=2, xaxt='n', yaxt='n', cex.lab=1.3,
log='y')

 #Estimates and confidence intervals
 points(c(1:4),c(1.1,1.32,1.14,1.36), pch=17, cex=1.5, col='blue')
 segments (c(1:4),c(0.93,1.11,0.94,1.15),c(1:4),c(1.3,1.58,1.37,1.61),
  col='blue', lwd=2)

 #Add x- and y-axes
 axis(1,c(1:4), labels= F)
 axis(2, at=seq(0.75,2, by=0.25), labels=seq(0.75,2, by=0.25), las=1)

 # Add x-axis labels
 labels - paste(Label, 1:4, sep =  )
 staxlab(side=1, at=1:4, labels, srt=45)

 -Original Message-
 From: Jim Lemon [mailto:j...@bitwrit.com.au]
 Sent: Tuesday, August 24, 2010 7:48 PM
 To: tesutton
 Cc: r-help@r-project.org
 Subject: Re: [R] Rotate x-axis label on log scale

 On 08/24/2010 11:44 AM, Tim Elwell-Sutton wrote:
 Hi

 I'd appreciate some help with plotting odds ratios. I want to  
 rotate the
 labels on the x-axis by 45 degrees.

 The usual way of doing this, using text - e.g. text(1,
 par('usr')[3]-2.25..)
 -  gives no result when the y-axis is a log scale.

 I guess this is because, as the par help says, for a logarithmic y- 
 axis:
 y-limits will be 10 ^ par(usr)[3:4]

 Hi Tim,
 If you know where to put the labels, try:

 library(plotrix)
 staxlab(1,at=...,labels=...,srt=45)

 Jim
-- 
David Winsemius, MD
West Hartford, CT

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


Re: [R] Odp: Finding pairs

2010-08-25 Thread Petr PIKAL
Hm

r-help-boun...@r-project.org napsal dne 25.08.2010 09:43:26:

 Dear Mr Petr Pikal
 
 I am extremely sorry for the manner I have raised the query. Actually 
that was
 my first post to this R forum and in fact even I was also bit confused 
while 
 drafting the query, for which I really owe sorry to all for consuming 
the 
 precious time. Perhaps I will try to redraft my query in a better way as 
follows. 
 
 I have two datasets A and B containing the names of branch offices 
of a 
 particular bank say XYZ plc bank. The XYZ bank has number of main branch 

 offices (say Parent) and some small branch offices falling under the 
purview 
 of these main branch offices (say Child).
 
 The datalist A and B consists of these main branch office names as 
well as
 small branch office names. B is subset of A and these branch names are 
coded. 
 Thus we have two datasets A and B as (again I am using only a
  portion of a large database just to have some idea)
 
 
 A B
 144  
   what is here in B? Empty space?, 
 145   
 146   
 147  144

How do you know that 144 from B relates to 147 in A? Is it according to 
its positions? I.e. 4th item in B belongs to 4.th item in A?

 148  145  
 
 149  147
 151  148
 
 
 
 Now the branch 144 appears in A as well as in B and in B it is mapped 
with 
 147. This means branch 147 comes under the purview of main branch 144. 
Again 
 147 is controlling the branch 149 (since 147 also has appeared in B and 
is 
 mapped with 149 of A).
 
 Similarly, branch 145 is controlling branch 148 which further controls 
 operations of bank branch 151 and like wise.

Well as you did not say anything about structure of your data
A-144:151
B-144:148
data.frame(A,B)
A   B
1 144  NA
2 145  NA
3 146  NA
4 147 144
5 148 145
6 149 146
7 150 147
8 151 148
DF-data.frame(A,B)
main-DF$A[is.na(DF$B)]
branch1-DF[!is.na(DF$B),]
selected.branch1-branch1$A[branch1$B%in%main]
branch2-branch1[!branch1$B%in%main,]
selected.branch2-branch2$A[branch2$B%in%selected.branch1]

and for cbinding your data which has uneven number of values see Jim 
Holtman's answer to this

How to cbind DF:s with differing number of rows?

Regards
Petr


 
 So in the end I need an output something like -
 
 Main Branch   Branch office1 Branch
  office2
 144 147 149
 145 148 151 
   
 146 NA
   NA   
 
...
 
..
 
  
 I understand again I am not able to put forward my query properly. But I 
must 
 thank all of you for giving a patient reading to my query and for 
reverting 
 back earlier. Thanks once again.
 
 With warmest regards
 
 Mike 
 
 
 --- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:
 
 From: Petr PIKAL petr.pi...@precheza.cz
 Subject: Odp: [R] Finding
  pairs
 To: Mike Rhodes mike_simpso...@yahoo.co.uk
 Cc: r-help@r-project.org
 Date: Wednesday, 25 August, 2010, 6:39
 
 Hi
 
 without other details it is probably impossible to give you any 
reasonable 
 advice. Do you have your data already in R? What is their form? Are they 

 in 2 columns in data frame? How did you get them paired?
 
 So without some more information probably nobody will invest his time as 

 it seems no trivial to me.
 
 Regards
 Petr
 
 r-help-boun...@r-project.org napsal dne 24.08.2010 20:28:42:
 
  
  
  
  
  Dear R Helpers,
  
  
  I am a newbie and recently got introduced to R. I have a large 
database 
  containing the names of bank branch offices along-with other details. 
I 
 am 
  into Operational
  Risk as envisaged by BASEL II Accord. 
  
  
  I am trying to express my problem and I am using only an indicative 
data 
 which
  comes in coded format.
  
  
  
  
  A (branch)  B (controlled by)
  
  
  144   
  145  
  146   
  147   144 
  148   145 
  149   
 147
  151   146  
   ..  ...
   
  ..  ...
  
  
  where 144's etc are branch codes in a given city and B is subset of A.
  
  
  
  
  If a branch code appearing in A also appears in B (which is paired 

 with 
  some otehr element of A e.g. 144 appearing in A, also appears in B 
and 
 is 
  paired with 147 of A and
  likewise), then that means 144 is controlling 
 
  operations of bank office 147. Again, 147 itself appears again 

[R] (no subject)

2010-08-25 Thread Mohammad Ali Vakili
Hi
I am using repeated meaturement data for my project and I want to use quantile 
regression with multilevel or panel data in R.
I dont find the basic version of software in R, so I have difficulty in using 
it.
I would also appreciate if anyone more proficient in R could help me how to run 
this.
best wishes 
M.A-Vakili

Department of Biostatistics 
Faculty of Medicine
Shiraz/Iran
P.O.Box: 71348-45794
E-Mail: vakil...@sums.ac.ir  or mavak...@yahoo.com



  
[[alternative HTML version deleted]]

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


Re: [R] Documenting S4 Methods

2010-08-25 Thread Martin Maechler
 Dario Strbenac d.strbe...@garvan.org.au
 on Wed, 25 Aug 2010 13:00:03 +1000 (EST) writes:

 I'm in the process of converting some S3 methods to S4 methods.
 I have this function :

 setGeneric(enrichmentCalc, function(rs, organism, seqLen, 
...){standardGeneric(enrichmentCalc)})

 setMethod(enrichmentCalc, c(GenomeDataList, BSgenome), function(rs, 
organism, seqLen, ...) {
 ...  ... ...
 })

 setMethod(enrichmentCalc, c(GenomeData, BSgenome), function(rs, 
organism, seqLen=NULL, do.warn=FALSE) {
 ...... ...
 })

 setMethod(enrichmentCalc, c(GRanges, BSgenome), function(rs, 
organism, seqLen=NULL) {
 ...... ...
 }

 and a part of my Rd file is :

 \name{enrichmentCalc}
 \docType{methods}
 \alias{enrichmentCalc,GenomeDataList,BSgenome-method}
 \alias{enrichmentCalc,GenomeData,BSgenome-method}
 \alias{enrichmentCalc,GRanges,BSgenome-method}
 ......   ...
 \usage{
 enrichmentCalc(rs, organism, seqLen, ...)
 enrichmentCalc(rs, organism, seqLen=NULL, do.warn=FALSE)
 enrichmentCalc(rs, organism, seqLen=NULL)
 }
 .........

 Can anyone suggest why I'm seeing this error :

 * checking for code/documentation mismatches ... WARNING
 Codoc mismatches from documentation object 'enrichmentCalc':
 enrichmentCalc
 Code: function(rs, organism, seqLen, ...)
 Docs: function(rs, organism, seqLen = NULL, do.warn = FALSE)
 Argument names in code not in docs:
 ...
 Argument names in docs not in code:
 do.warn
 Mismatches in argument names:
 Position: 4 Code: ... Docs: do.warn
 Mismatches in argument default values:
 Name: 'seqLen' Code:  Docs: NULL
 enrichmentCalc
 Code: function(rs, organism, seqLen, ...)
 Docs: function(rs, organism, seqLen = NULL)
 Argument names in code not in docs:
 ...
 Mismatches in argument default values:
 Name: 'seqLen' Code:  Docs: NULL

 * checking Rd \usage sections ... WARNING
 Objects in \usage without \alias in documentation object 'enrichmentCalc':
 enrichmentCalc

 Also, what is the difference between

 ...  ...  ...
 \docType{methods}
 ...  ...  ...
 \alias{methodName,class-method}
 ...  ...  ...
 \usage{methodName(arg1)}
 ...  ...  ...

 and

 ...  ...  ...
 \alias{methodName,class-method}
 ...  ...  ...
 \usage
 {
 \S4method{methodName}{class}(arg1)
 }
 ...  ...  ...

the last one is the one you should use.

I don't think that you can easily use \usage{..} without
\S4method{.}
and not get a warning (/error) from R CMD check.


 I've seen both ways used for S4 methods and don't know what is the 
underlying difference.

BTW: You should probably reread (parts of) the Writing R Extensions manual;
 It's a good exercise ... even for me as an R core team member.


 I haven't been able to find any good tutorials for the new S4 
architecture (written post 2006), so I'm not sure where to start with S4.

In R,  ?Methods  - References  (s r in ESS)
lists John Chambers book (with extra hints).

 Chambers, John M. (2008) _Software for Data Analysis: Programming
 with R_ Springer.  (For the R version: see section 10.6 for method
 selection and section 10.5 for generic functions).

Many people recommend learning from good examples.
As the Bioconductor project has made heavy use of S4 classes and
methods, many of their packages are good examples to follow.

To find all - locally installed - packages which depend directly
on methods and hence are probably using S4 methods / classes,
I do

 library(tools)
 str(meth.pkgs.direct - dependsOnPkgs(methods, recursive=FALSE))
 chr [1:313] BufferedMatrix externalVector gdistancetest Matrix ...
 sort(meth.pkgs.direct)
  [1] AcceptanceSamplingaccuracy  adegenet 
  [4] adephylo  amer  aod  
  [7] apcluster apsrtable archetypes   
 [10] ArDec arulesarulesNBMiner
 [13] arulesSequences   asuR  automap  
 [16] aws   aylmerbackfitRichards  
 [19] backtest  bbmle bcp  
 [22] bdsmatrix biclust   bifactorial  
 [25] biganalytics  biglm bigmemory
 [28] bigtabulate   bild  BLCOP
 [31] Brobdingnag   bsBufferedMatrix   
 [34] calib catnetcelsius  
 [37] ChainLadder   classGraphclue 
 [40] 

Re: [R] how to plot y-axis on the right of x-axis

2010-08-25 Thread Jim Lemon

On 08/25/2010 09:12 AM, elaine kuo wrote:

Dear List,

I have a richness data distributing across 20 N to 20 S latitude. (120 E-140
E longitude).


I would like to draw the richness in the north hemisphere and a regression
line in the plot
(x-axis: latitude, y-axis: richness in the north hemisphere).
The above demand is done using plot.

Then, south hemisphere richness and regression are required to be generated
using
the same y-axis above but an x-axis on the left side of the y-axis.
(The higher latitude in the south hemisphere, the left it would move)

Please kindly share how to design the south plot and regression line for
richness.
Also, please advise if any more info is in need.


Hi Elaine,
Changing the position of the axes is easily done by changing the side 
and pos arguments of the axis function. If you want to move the 
y-axis to the right (or as you state it, the x axis to the left):


# y axis on the left
plot(1:5,axes=FALSE)
axis(1)
axis(2)
# add a y axis one third of the way to the right
xylim-par(usr)
axis(2,pos=xylim[1]+diff(xylim[1:2])/3)
# add another y axis two thirds of the way
axis(4,pos=xylim[2]-diff(xylim[1:2])/3)
# add one more on the right
axis(4)

You can move the x axis up and down using the same tricks.

Jim

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


Re: [R] percentage sign in expression

2010-08-25 Thread Gavin Simpson
On Wed, 2010-08-25 at 09:32 +0100, e-letter wrote:
 On 24/08/2010, David Winsemius dwinsem...@comcast.net wrote:
 
  On Aug 24, 2010, at 9:37 AM, e-letter wrote:
 
  Readers,
 
  According to the documentation for the function 'plotmath' there is no
  apparent possibility to add the percent sign (%) to a plot function,
 
  Where did you see an assertion made???
 
 Within R I entered the command:
 
 ?plotmath

Surely you don't expect the R help pages to document everything that
*isn't* possible with a given function? They'd be infinitely long! ;-)

What you meant to say was, ?plotmath doesn't show me how to add a %
symbol to a plot. How do I do it?

 Also accessed using:
 
 help.start(browser=opera)
 
 Navigating the web browser page:
 
 packages
 packages in /usr/lib/R/library
 grdevices
 plotmath
 
 In the list headed 'syntax' and 'meaning' within the section 'details'.
 
snip /
TestChars - function(sign=1, font=1, ...)
  + {
  +if(font == 5) { sign - 1; r - c(32:126, 160:254)
  +} else if (l10n_info()$MBCS) r - 32:126 else r - 32:255
  +if (sign == -1) r - c(32:126, 160:255)
  +par(pty=s)
  +plot(c(-1,16), c(-1,16), type=n, xlab=, ylab=,
  + xaxs=i, yaxs=i)
  +grid(17, 17, lty=1)
  +for(i in r) try(points(i%%16, i%/%16, pch=sign*i, font=font,...))
  + }
TestChars(font=5)
 
  Notice that the % sign is three characters to the right (i.e.
  higher) of the forall symbol that is produced by the example code
 
 I can't see 'forall' in the code above.

You need to run it in R. Then you will see a plot of the glyphs
available in the symbol font. 'forall' is a mathematical symbol:

http://en.wikipedia.org/wiki/Table_of_mathematical_symbols

  they offer. (The numbering proceeds from bottom up which confused me
  at first.)
 
 What numbering?

The numbering of the glyphs so you can use the number to draw the symbol
you want. They are on the plot. You did run the code provided by David?

Anyway, is there a problem with this:

plot(1:10, xlab = expression(foo == 15*%))

That way you don't need to worry about finding the right glyph.

HTH

G

 The documentation makes reference to the command:
 
 demo(plotmath)
 
 I applied this command and could not see an instruction to produce the
 percent (%) symbol.
 
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,  [f] +44 (0)20 7679 0565
 Pearson Building, [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London  [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT. [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

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


Re: [R] Cran-packages-ProfessoR- how to Automatically email a file to an address using the perl program.

2010-08-25 Thread veepsirtt

Hi Jonathan lees

How to use this code and after installing ProfessR  package

autoemail(eadd, sfile, hnote = Exam Results) 

thanks
veepsirtt
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Cran-packages-ProfessoR-how-to-Automatically-email-a-file-to-an-address-using-the-perl-program-tp2335184p2337924.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] help

2010-08-25 Thread Ingmar Visser
Martin,

On Wed, Aug 25, 2010 at 10:17 AM, Loos, Martin martin.l...@eawag.ch wrote:

 Good afternoon!

 It may seem trivial to some/most of You, but I found it difficult to
 properly include a C++-based .dll into a package that I want to build for
 usage in R. I read through the Writing R extensions...  R administration
 ... instructions, but it seems I did not grasp the bigger picture.

 The way I figured out to use the .dll in my package finally worked for
 using that package from the R console, but it gives multiple error and
 warning messages when running the RCMD check pkg command in the shell.

 What I did: I created the package skeleton via package.skeleton() including
 a function that makes both (1) a
 library.dynam() call to a .dll and (2) a .C() call including the name of
 the function in that DLL. Next, I run RCMD INSTALL pkg to install the
 package and then created a /lib folder in the installed package containing
 the .dll file.


As you say below, the way to go here is to add a /src folder in the package
skeleton which contains your C/C++ code and
header files and then run R CMD INSTALL, which will compile your code into a
.dll file and place it in the appropriate directory,
ie the /lib directory. It should not be necessary to do this manually.

hth, Ingmar

PS: It's hard to tell what is going wrong without the actual error messages.




 I tried alternatives, e.g., (a) containing a /src folder in may package
 skeleton with the .dll in it and/or (b) including zzz.R with a .First.lib
 function including the library call  but nothing worked out when running
 the package check.

 help! thank You very much, MArtin

[[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] Rotate x-axis label on log scale

2010-08-25 Thread Jim Lemon

On 08/25/2010 01:37 PM, Tim Elwell-Sutton wrote:

Hi Jim
Thanks for this. The staxlab function seems very useful. Unfortunately, the
rotation option doesn't seem to work for me when the y-axis is on a log
scale. It will stagger the labels but not rotate them. There's no error
message. On a linear axis the rotation works nicely. Any ideas?
The example below works if you omit log='y' or srt=45


Hi Tim,
Thanks for letting me know about this problem. I never did get around to 
making staxlab work with log axes, but I think this new version:


staxlab-function(side=1,at,labels,nlines=2,top.line=0.5,
 line.spacing=0.8,srt=NA,ticklen=0.03,...) {

 if(missing(labels)) labels-at
 nlabels-length(labels)
 if(missing(at)) at-1:nlabels
 if(is.na(srt)) {
  linepos-rep(top.line,nlines)
  for(i in 2:nlines) linepos[i]-linepos[i-1]+line.spacing
  linepos-rep(linepos,ceiling(nlabels/nlines))[1:nlabels]
  axis(side=side,at=at,labels=rep(,nlabels))
  mtext(text=labels,side=side,line=linepos,at=at,...)
 }
 else {
  xylim-par(usr)
  if(side == 1) {
   xpos-at
   if(par(ylog)) ypos-10^(xylim[3]-ticklen*(xylim[4]-xylim[3]))
   else ypos-xylim[3]-ticklen*(xylim[4]-xylim[3])
  }
  else {
   ypos-at
   if(par(xlog)) xpos-10^(xylim[1]-ticklen*(xylim[2]-xylim[1]))
   else xpos-xylim[1]-ticklen*(xylim[2]-xylim[1])
  }
  par(xpd=TRUE)
  text(xpos,ypos,labels,srt=srt,adj=1,...)
  par(xpd=FALSE)
 }
}

will do what you want.

Jim

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


Re: [R] Plot bar lines like excel

2010-08-25 Thread Dennis Murphy
Hi:

Here's a ggplot2 version of your graph. Since you didn't include the dates,
I had to make them up, and it may have some consequence with respect to your
times because yours don't appear to be exactly equally spaced.

I created two data frames - an 'original' one that contains both series in
separate columns, and another that 'melts' the original data frame by
stacking the two series horizontally with a created factor that holds the
variable names. Notice that I assigned the names of the two series to
Series4 and Series6 in the original data frame to correspond with the names
on the plot.

library(ggplot2)
dd - data.frame(tm = seq(as.Date('2002-1-1'), by = 'month', length = 50),
 Series4 = x1, Series6 = x2)
dm - melt(dd, id = 'tm')
dd - transform(dd, ylo = pmin(Series4, Series6), yup = pmax(Series4,
Series6))

The melted data, dm, is useful for plotting the two series and associating a
color aesthetic to the different series; we use it to generate the initial
part of the plot, the two series themselves. 'variable' is a factor variable
whose levels are the variable names (which map to the legend) and 'values'
is a numeric variable containing the values of the two series.

g - ggplot(data = dm, aes(x = tm))
g + geom_line(aes(y = value, colour = variable), size = 1.2) +
scale_colour_manual(, values = c('Series4' = 'yellow',
   'Series6' = 'red')) + theme_bw() + ylim(4, 12) +
geom_linerange(data = dd, aes(x = tm, ymin = ylo, ymax = yup), size =
0.8) +
opts(panel.grid.major = theme_blank()) +
opts(panel.grid.minor = theme_blank()) +
geom_hline(yintercept = c(4, 6, 8, 10, 12), size = 1, alpha = 0.1) +
scale_x_date(major = '2 months', format = '%b-%Y') +
opts(axis.text.x = theme_text(angle = 90, hjust = 1, vjust = 0.5, size =
9))

scale_colour_manual() allows specification of one's own colors rather than
ggplot2's defaults. The initial pair of double quotes gets rid of the legend
label - if you leave it off, it will use the header 'variable'. You can
always substitute in another legend header if you wish.

The vertical lines connecting the series at each date is where the work is
required. ylo and yup in data frame dd represent the min and max values of
the series at each date, respectively, produced by the pmin() and pmax()
functions. ggplot2's geom_linerange() needs these variables to represent the
lower and upper limits of the line segments at each x value (date), which is
why the data frame used in geom_linerange() is dd rather than dm. The size
parameter in geom_line() and geom_linerange() corresponds to line thickness.


The first two opts() statements get rid of the default grid lines;
geom_hline() provides the substitutes.
The scale_x_date() code basically labels every other date with some
additional options to pretty up the appearance of the labels.

HTH,
Dennis

On Tue, Aug 24, 2010 at 2:40 PM, abotaha yaseen0...@gmail.com wrote:


 Thanks a lot for the nice explanation. however, you mention that it will be
 nicer if there is data. yes you are right it be very kind of you if you
 include a simple example. where my data are:

 x1-c(11.5,9.38,9.3,9.8,9,9.06,9.42,8.8,9.05,8.14,8.2,7.59,6.92,6.47,7.12,

 7.47,6.81,8.41,9.64,9.62,9.2,8.92,8,7.6,7.6,7.19,6.94,6.91,7,7.25,6.6,7.18,7.78,

 7.37,7.29,6.71,7.05,6.69,6.05,6.38,6.18,7.67,7.34,7.13,6.95,6.8,6.45,6.81,6.27,6.35)


 x2-c(11.5,8.959535,9.728104,9.609262,8.755494,9.221545,8.244458,7.63944,7.92052,7.690449,7.853247,

 7.239616,6.609007,5.92946,6.47822,5.906936,6.966004,8.630517,9.506582,9.57479,9.236638,9.581875,8.838992,8.368731,8.608884,7.998023,7.918123,7.832322,7.930177,7.479222,6.866978,7.454062,8.206185,7.344037,7.059774,6.547646,7.005803,6.623987,5.992691,6.313481,6.194613,6.224266,7.084932,6.568976,6.43866,5.70081,6.7593,6.6929,6.46806,5.964816)

 regards,



 --
 View this message in context:
 http://r.789695.n4.nabble.com/Plot-bar-lines-like-excel-tp2337089p2337394.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.


[[alternative HTML version deleted]]

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


Re: [R] Plot bar lines like excel

2010-08-25 Thread abotaha

Woow, it is amazing, 
thank you very much. 
yes i forget to attach the dates, however, the dates in my case is every 16
days. 
so how i can use 16 day interval instead of month in by option.

cheers, 


-- 
View this message in context: 
http://r.789695.n4.nabble.com/Plot-bar-lines-like-excel-tp2337089p2338028.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] Odp: Finding pairs

2010-08-25 Thread Mike Rhodes
Dear Mr Petr PIKAL
After reading the R code provided by you, I realized that I would have never 
figured out how this could have been done. I am going to re-read again and 
again your code to understand the logic and the commands you have provided.
Thanks again from the heart for your kind advice.
Regards
Mike

--- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:

From: Petr PIKAL petr.pi...@precheza.cz
Subject: Re: [R] Odp:  Finding pairs
To: Mike Rhodes mike_simpso...@yahoo.co.uk
Cc: r-help@r-project.org
Date: Wednesday, 25 August, 2010, 9:01

Hm

r-help-boun...@r-project.org napsal dne 25.08.2010 09:43:26:

 Dear Mr Petr Pikal
 
 I am extremely sorry for the manner I have raised the query. Actually 
that was
 my first post to this R forum and in fact even I was also bit confused 
while 
 drafting the query, for which I really owe sorry to all for consuming 
the 
 precious time. Perhaps I will try to redraft my query in a better way as 
follows. 
 
 I have two datasets A and B containing the names of branch offices 
of a 
 particular bank say XYZ plc bank. The XYZ bank has number of main branch 

 offices (say Parent) and some small branch offices falling under the 
purview 
 of these main branch offices (say Child).
 
 The datalist A and B consists of these main branch office names as 
well as
 small branch office names. B is subset of A and these branch names are 
coded. 
 Thus we have two datasets A and B as (again I am using only a
  portion of a large database just to have some idea)
 
 
 A                         B
 144                      
                       what is here in B? Empty space?, 
 145                       
 146                       
 147                  144                        

How do you know that 144 from B relates to 147 in A? Is it according to 
its positions? I.e. 4th item in B belongs to 4.th item in A?

 148                  145  
 
 149                  147
 151                  148
 
 
 
 Now the branch 144 appears in A as well as in B and in B it is mapped 
with 
 147. This means branch 147 comes under the purview of main branch 144. 
Again 
 147 is controlling the branch 149 (since 147 also has appeared in B and 
is 
 mapped with 149 of A).
 
 Similarly, branch 145 is controlling branch 148 which further controls 
 operations of bank branch 151 and like wise.

Well as you did not say anything about structure of your data
A-144:151
B-144:148
data.frame(A,B)
    A   B
1 144  NA
2 145  NA
3 146  NA
4 147 144
5 148 145
6 149 146
7 150 147
8 151 148
DF-data.frame(A,B)
main-DF$A[is.na(DF$B)]
branch1-DF[!is.na(DF$B),]
selected.branch1-branch1$A[branch1$B%in%main]
branch2-branch1[!branch1$B%in%main,]
selected.branch2-branch2$A[branch2$B%in%selected.branch1]

and for cbinding your data which has uneven number of values see Jim 
Holtman's answer to this

How to cbind DF:s with differing number of rows?

Regards
Petr


 
 So in the end I need an output something like -
 
 Main Branch           Branch office1                 Branch
  office2
 144                             147                                 149
 145                             148                                 151 
   
 146                             NA
                                   NA               
 
...
 
..
 
  
 I understand again I am not able to put forward my query properly. But I 
must 
 thank all of you for giving a patient reading to my query and for 
reverting 
 back earlier. Thanks once again.
 
 With warmest regards
 
 Mike 
 
 
 --- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:
 
 From: Petr PIKAL petr.pi...@precheza.cz
 Subject: Odp: [R] Finding
  pairs
 To: Mike Rhodes mike_simpso...@yahoo.co.uk
 Cc: r-help@r-project.org
 Date: Wednesday, 25 August, 2010, 6:39
 
 Hi
 
 without other details it is probably impossible to give you any 
reasonable 
 advice. Do you have your data already in R? What is their form? Are they 

 in 2 columns in data frame? How did you get them paired?
 
 So without some more information probably nobody will invest his time as 

 it seems no trivial to me.
 
 Regards
 Petr
 
 r-help-boun...@r-project.org napsal dne 24.08.2010 20:28:42:
 
  
  
  
  
  Dear R Helpers,
  
  
  I am a newbie and recently got introduced to R. I have a large 
database 
  containing the names of bank branch offices along-with other details. 
I 
 am 
  into Operational
  Risk as envisaged by BASEL II Accord. 
  
  
  I am trying to express my problem and I am using only an indicative 
data 
 which
  comes in coded format.
  
  
  
  
  A (branch)                      B (controlled by)
  
  
  144                   
  145                      
  146                   
  147                                       144 
  148                                       145 
  149       
                                

[R] Merging two data set in R,

2010-08-25 Thread Mangalani Peter Makananisa
Dear R Gurus,

 

I am currently working on the two dataset ( A and B), they both have the
same fields:ID , REGION, OFFICE, CSTART, CEND, NCYCLE, STATUS and
CB.

I want to merge the two data set by ID. The problem I have is that the
in data A, the ID's are unique. However in the data set B, the ID's are
not unique, thus some repeat themselves.

 

How do I the merge or retrieve the common ones? 

Please advise.

 

Kind Regards

 

Peter

 

South Africa

+27 12 422 7357

+27 82 456 4669

 

 

 

 

 

 

 


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

[[alternative HTML version deleted]]

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


[R] keys in score.items in package psych

2010-08-25 Thread Frans Marcelissen

Hi Ruru's,
score.items in package psych is very handy for scoring test items. It has the 
structure score.items(keys,items).
For instance: 
score.items(c(1,1,1),data.frame(a1=rep(1,5),a2=rep(1,5),a3=rep(1,5)))$scores
correctly gives 1 on each case.
But if key  -1,0,1  the following happens
score.items(c(0.1,0.1,0.1),data.frame(a1=rep(1,5),a2=rep(1,5),a3=rep(1,5)))$scores
gives 10 on each case. It appears that the inverse of key is the relative 
weight of that item.
My questions: is this correct? Is this an undocumented feature? Can i trust on 
it?
Frans
  
[[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] Cran-packages-ProfessoR- how to Automatically email a file to an address using the perl program.

2010-08-25 Thread Jonathan Lees-2

Yes -
please look at the function IDandEM.R
(This stand for IDentification and EMail)

in there it has a complete function that uses
the autoemail function - I use this all the time.

There is a logical  option (SEND)
on whether to actually mail the
file to the recipient.

The IDandEM function compares two lists
(one from the exam the other from
the list that stores the email addresses) and matches these before mailing
the scores to the students.  I do this to help prevent accidentally
mailing students the wrong grades.
In the U.S. that would be a very serious breech of protocol.

I usually do not run the autoemail directly because of this.



-- 
Jonathan M. Lees
Professor
THE UNIVERSITY OF NORTH CAROLINA AT CHAPEL HILL
UNC-CH Geological Sciences
Mitchell Hall
104 South Road
CB 3315
Chapel Hill, NC  27599-3315
TEL: (919) 962-1562
FAX: (919) 966-4519
jonathan.l...@unc.edu
http://www.unc.edu/~leesj



-- 
View this message in context: 
http://r.789695.n4.nabble.com/Cran-packages-ProfessoR-how-to-Automatically-email-a-file-to-an-address-using-the-perl-program-tp2335184p2337984.html
Sent from the R help mailing list archive at Nabble.com.

[[alternative HTML version deleted]]

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


[R] using an objects contents in a text string

2010-08-25 Thread josquint

Hi,  I am (trying) to write a script that will execute a set of functions and
then to write some of the output data frames to file. I am wanting to be
able to just change a single object and have this them populate through the
functions etc. I have no trouble making this work when the function is
calling for an argument that meets this criteria

e.g.

input.variable.name - DDDHHD

list - function(x, input.variable = input.variable.name, ...)

where DDDHHD will be then used as the input.variable in the function BUT

I can not seem to get the text string used as part of the name of the output
file name in the write.data.frame function.

e.g.

write(x, file = input.variable.name _data.txt,  ...)

to create an output file with the name DDDHHD_data.txt

I know this syntax is incorrect but could someone guide me as to how to
acieve this (I hope it is clear??).

Thanks in advance.


-- 
View this message in context: 
http://r.789695.n4.nabble.com/using-an-objects-contents-in-a-text-string-tp2338061p2338061.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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread skan

 down vote  favorite


Hello

I have a zoo series. It lasts 10 years and its frequency is 15min.

I'd like to get a new zoo series (or vector) with the same number of
elements, whith each element equal to the first element of the day. That's,
The first element everyday is repeated throughout the wole day.

This is not same as aggregate(originalseries,as.Date,head,1) because this
gives a vector with just one element for each day.

cheers

-- 
View this message in context: 
http://r.789695.n4.nabble.com/Repeat-the-first-day-data-through-all-the-day-Zoo-tp2338069p2338069.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] lmer() causes segfault

2010-08-25 Thread Dennis Murphy
Hi:

Let's start with the data:

 str(test.data)
'data.frame':   100 obs. of  4 variables:
 $ StudentID: num  17370 17370 17370 17370 17379 ...
 $ GroupID  : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Time : num  1 2 3 4 1 2 3 4 1 2 ...
 $ Score: num  76.8 81.8 89.8 92.8 75.9 ...

Both StudentID and GroupID are numeric; in the model, they would be treated
as continuous covariates rather than factors, so we need to convert:

test.data$StudentID - factor(test.data$StudentID)
test.data$GroupID - factor(test.data$GroupID)

Secondly, I believe there are some flaws in your model. After converting
your variables to factors, I ran

library(lme4)
mlmoded1.lmer - lmer(Score ~ Time + (Time | GroupID/StudentID), data =
test.data)

You have two groups, so they should be treated as a fixed effect - more
specifically, as a fixed blocking factor. The StudentIDs are certainly
nested within GroupID, and Time is measured on each StudentID, so it is a
repeated measures factor. The output of this model is
 mlmoded1.lmer
Linear mixed model fit by REML
Formula: Score ~ Time + (Time | GroupID/StudentID)
   Data: test.data
   AIC   BIC logLik deviance REMLdev
 393.1 416.5 -187.5376.9   375.1
Random effects:
 GroupsNameVariance  Std.Dev. Corr
 StudentID:GroupID (Intercept)  0.504131 0.71002
   Time 0.083406 0.28880  1.000
 GroupID   (Intercept) 12.809567 3.57905
   Time 3.897041 1.97409  -1.000
 Residual   1.444532 1.20189
Number of obs: 100, groups: StudentID:GroupID, 25; GroupID, 2

Fixed effects:
Estimate Std. Error t value
(Intercept)   72.803  2.552  28.530
Time   4.474  1.401   3.193

Correlation of Fixed Effects:
 (Intr)
Time -0.994

The high correlations among the random effects and then among the fixed
effects suggests that the model specification may be a bit off.

The above model fits random slopes to GroupIDs and StudentIDs, along with
random intercepts, but GroupID is a between-subject effect and should be at
the top level. Time is a within-subject effect and StudentIDs are the
observational units. I modified the model to provide fixed effects for
GroupIDs, scalar random effects for StudentIDs and random slopes for
StudentIDs.

 mod3 - lmer(Score ~ 1 + GroupID + Time + (1 | StudentID) +
+  (0 + Time | StudentID), data = test.data)
 mod3
Linear mixed model fit by REML
Formula: Score ~ 1 + GroupID + Time + (1 | StudentID) + (0 + Time |
StudentID)
   Data: test.data
   AIC   BIC logLik deviance REMLdev
 430.9 446.5 -209.4418.4   418.9
Random effects:
 GroupsNameVariance   Std.Dev.
 StudentID (Intercept) 4.2186e-13 6.4951e-07
 StudentID Time1.8380e+00 1.3557e+00
 Residual  1.6301e+00 1.2768e+00
Number of obs: 100, groups: StudentID, 25

Fixed effects:
Estimate Std. Error t value
(Intercept)  70.7705 0.4204  168.33
GroupID2  4.0248 0.58546.88
Time  4.5292 0.2942   15.39

Correlation of Fixed Effects:
 (Intr) GrpID2
GroupID2 -0.668
Time -0.264  0.000

I didn't check the quality of the fit, but on the surface it seems to be
more stable, FWIW. Perhaps one could also add a term (GroupID | StudentID),
but I don't know offhand if that would make any sense. Another issue to
consider is whether to fit by REML or ML, but that is secondary to getting
the form of the model equation right. I don't claim this as a final model,
but rather a 're-starting point'. It may well be in need of improvement, so
comments are welcome.

The confusion between subjects nested in time or vice versa has occurred
several times this week with respect to repeated measures/longitudinal
models using lmer(), so perhaps it merits a comment: subjects/experimental
units are NOT nested in time. Measurements taken on an individual at several
time points *entails* that time be nested within subject. Just saying...

This discussion may be better continued on the R-sig-mixed list, so I've
cc-ed to that group as well.

HTH,
Dennis

On Wed, Aug 25, 2010 at 1:27 AM, Bertolt Meyer bme...@sozpsy.uzh.ch wrote:

 Ben Bolker bbolker at gmail.com writes:

  Bertolt Meyer bmeyer at sozpsy.uzh.ch writes:


 Hello lmer() - users,

 A call to the lmer() function causes my installation of R (2.11.1 on
 Mac OS X 10.5.8) to crash and I am trying to figure out the problem.


  [snip snip]

  detach(package:nlme)
 library(lme4)

 mod1 - lmer(performance ~ time + (time | GroupID/StudentNumber), data
 = dataset.long, na.action = na.omit)

 However, this call results in a segfault:

 *** caught segfault ***
 address 0x154c3000, cause 'memory not mapped'

 and a lengthy traceback. I can reproduce this error. It also occurs
 when I don't load nlme before lme4. Can someone tell me what I am
 doing wrong? Any help is greatly appreciated.


  This may well be a bug in lmer.  There have been a number of
 fussy computational issues with the lme4 package on the Mac 

[R] Surprising behaviour survey-package with missing values

2010-08-25 Thread Jan van der Laan
Dear list,

I got some surprising results when using the svytotal routine from the
survey package with data containing missing values.

Some example code demonstrating the behaviour is included below.

I have a stratified sampling design where I want to estimate the total
income. In some strata some of the incomes are missing. I want to
ignore these missing incomes. I would have expected that
svytotal(~income, design=mydesign, na.rm=TRUE) would do the trick.
However, when calculating the estimates 'by hand' the estimates were
different from those obtained from svytotal. The estimated mean
incomes do agree with each other. It seems that using the na.rm option
with svytotal is the same as replacing the missing values with zero's,
which is not what I would have expected, especially since this
behaviour seems to differ from that of svymean. Is there a reason for
this behaviour?

I can of course remove the missing values myself before creating the
survey object. However, with many different variables with different
missing values, this is not very practical. Is there an easy way to
get the behaviour I want?

Thanks for your help.

With regards,

Jan van der Laan


=== EXAMPLE ===

library(survey)
library(plyr)

# generate some data
data - data.frame(
  id = 1:20,
  stratum = rep(c(a, b), each=10),
  income = rnorm(20, 100),
  n = rep(c(100, 200), each=10)
  )
data$income[5] - NA

# calculate mean and total income for every stratum using survey package
des - svydesign(ids=~id, strata=~stratum, data=data, fpc=~n)
svyby(~income, by=~stratum, FUN=svytotal, design=des, na.rm=TRUE)
mn  - svyby(~income, by=~stratum, FUN=svymean, design=des, na.rm=TRUE)
mn
n   - svyby(~n, by=~stratum, FUN=svymean, design=des)

# total does not equal mean times number of persons in stratum
mn[2] * n[2]

# calculate mean and total income 'by hand'. This does not give the same total
# as svytotal, but it does give the same mean
ddply(data, .(stratum), function(d) {
data.frame(
  mean = mean(d$income, na.rm=TRUE),
  n = mean(d$n),
  total = mean(d$income, na.rm=TRUE) * mean(d$n)
)
  })

# when we set income to 0 for missing cases and repeat the previous estimation
# we get the same answer as svytotal (but not svymean)
data2 - data
data2$income[is.na(data$income )] - 0
ddply(data2, .(stratum), function(d) {
data.frame(
  mean = mean(d$income, na.rm=TRUE),
  n = mean(d$n),
  total = mean(d$income, na.rm=TRUE) * mean(d$n)
)
  })

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

2010-08-25 Thread Sarah Goslee
First you need to clarify what you'd like to happen when the ID in B
is not unique. What do you want the resulting dataframe to look
like?

Some possible answers involve using different options for merge() or
using unique() to remove duplicates from B before merging. But
at least to me, merge or retrieve the common ones isn't clear
enough to be able to say which.

Sarah

On Wed, Aug 25, 2010 at 5:35 AM, Mangalani Peter Makananisa
pmakanan...@sars.gov.za wrote:
 Dear R Gurus,



 I am currently working on the two dataset ( A and B), they both have the
 same fields:    ID , REGION, OFFICE, CSTART, CEND, NCYCLE, STATUS and
 CB.

 I want to merge the two data set by ID. The problem I have is that the
 in data A, the ID's are unique. However in the data set B, the ID's are
 not unique, thus some repeat themselves.



 How do I the merge or retrieve the common ones?

 Please advise.



 Kind Regards



 Peter


-- 
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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread Gabor Grothendieck
On Wed, Aug 25, 2010 at 7:43 AM, skan juanp...@gmail.com wrote:
 I have a zoo series. It lasts 10 years and its frequency is 15min.

 I'd like to get a new zoo series (or vector) with the same number of
 elements, whith each element equal to the first element of the day. That's,
 The first element everyday is repeated throughout the wole day.

 This is not same as aggregate(originalseries,as.Date,head,1) because this
 gives a vector with just one element for each day.

Try ave:

 library(zoo)
 library(chron)
 zz - z - zoo(1:100, chron(0:9/5))
 zz[] - ave(coredata(z), as.Date(time(z)), FUN = function(x) head(x, 1))
 cbind(z, zz)
 z zz
(01/01/70 00:00:00)  1  1
(01/01/70 04:48:00)  2  1
(01/01/70 09:36:00)  3  1
(01/01/70 14:24:00)  4  1
(01/01/70 19:12:00)  5  1
(01/02/70 00:00:00)  6  6
(01/02/70 04:48:00)  7  6
(01/02/70 09:36:00)  8  6
(01/02/70 14:24:00)  9  6
(01/02/70 19:12:00) 10  6

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] multi day R run crashed - is there a log

2010-08-25 Thread Sarah Goslee
Oh dear. No there isn't. I do a lot of very long runs, and have learned to write
out intermediate steps. Depending on what you are doing, saving the RData
file periodically may be appropriate, or writing out a csv file of
results so far,
or even just printing something to the output file (if you are using batch).

There may well be something more elegant, but these solutions have all worked
well for me in various situations.

Sarah

On Wed, Aug 25, 2010 at 3:46 AM, Martin Tomko martin.to...@geo.uzh.ch wrote:
 Dear all,
 I am using an R 2.10 installation on a Windows 203 server that I have no
 control over. After a multi-day run I found that it was terminated/crashed.
 Is there any log kept by R where I could see whether something/what
 happened? The same process has been run beofre on a smaller dataset (also at
 least a day of computing) without problems.

 Thanks
 Martin

-- 
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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread Gabor Grothendieck
On Wed, Aug 25, 2010 at 7:56 AM, Gabor Grothendieck
ggrothendi...@gmail.com wrote:
 On Wed, Aug 25, 2010 at 7:43 AM, skan juanp...@gmail.com wrote:
 I have a zoo series. It lasts 10 years and its frequency is 15min.

 I'd like to get a new zoo series (or vector) with the same number of
 elements, whith each element equal to the first element of the day. That's,
 The first element everyday is repeated throughout the wole day.

 This is not same as aggregate(originalseries,as.Date,head,1) because this
 gives a vector with just one element for each day.

 Try ave:

 library(zoo)
 library(chron)
 zz - z - zoo(1:100, chron(0:9/5))

That should have been 10, not 100; however, it ignored 11:100 so the
answer is the same.

 zz[] - ave(coredata(z), as.Date(time(z)), FUN = function(x) head(x, 1))
 cbind(z, zz)
                     z zz
 (01/01/70 00:00:00)  1  1
 (01/01/70 04:48:00)  2  1
 (01/01/70 09:36:00)  3  1
 (01/01/70 14:24:00)  4  1
 (01/01/70 19:12:00)  5  1
 (01/02/70 00:00:00)  6  6
 (01/02/70 04:48:00)  7  6
 (01/02/70 09:36:00)  8  6
 (01/02/70 14:24:00)  9  6
 (01/02/70 19:12:00) 10  6


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

2010-08-25 Thread Sarah Goslee
Hello,

http://www.rseek.org is extremely useful if you know what you want to do but
not how to do it. Going there and putting quantile regression with multilevel
data into the search box returns many references to functions and packages
that may be of use.

This list can't do much more for you without a clearer statement of your
problem, as requested by the posting guide linked at the bottom of this
message.

Also, for future reference: I suspect many people on this list just delete
messages with no subject. You will receive faster and more informed
advice if you use a meaningful subject line.

On Wed, Aug 25, 2010 at 5:32 AM, Mohammad Ali Vakili mavak...@yahoo.com wrote:
 Hi
 I am using repeated meaturement data for my project and I want to use quantile
 regression with multilevel or panel data in R.
 I dont find the basic version of software in R, so I have difficulty in using
 it.
 I would also appreciate if anyone more proficient in R could help me how to 
 run
 this.
 best wishes
 M.A-Vakili



-- 
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] using an objects contents in a text string

2010-08-25 Thread josquint

ALSO    I have had a play with cat() but have also not got this to
work


e.g. 

write(x, file = cat(input.variable.name  , file = , sep = _data.txt, ),
...)

but this does not seem to work and I'm sure it is not the correct use of
cat()

Thanks.
-- 
View this message in context: 
http://r.789695.n4.nabble.com/using-an-objects-contents-in-a-text-string-tp2338061p2338106.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] Merging two data set in R,

2010-08-25 Thread Sarah Goslee
What do you want to happen when there are duplicates?

A:
ID X
1  a
2  b
3  c

B:
ID Y
1  x
2  y
2  z

What happens to ID 1? 2? 3? in your desired output?

The all.x and all.y options might be of use.

Sarah

On Wed, Aug 25, 2010 at 8:00 AM, Mangalani Peter Makananisa
pmakanan...@sars.gov.za wrote:
 I want to merge data set A and B, by merge(A,B, by = ID), however I am 
 getting  error massages, because the some ID's in A repeat themselves several 
 time in data set B. Even if the ID's in B repeat themselves I want to be able 
 to merge the two dataset and retrieve the intersection.

 Please help.

 -Original Message-
 From: Sarah Goslee [mailto:sarah.gos...@gmail.com]
 Sent: 25 August 2010 01:52 PM
 To: Mangalani Peter Makananisa
 Cc: r-help@r-project.org
 Subject: Re: [R] Merging two data set in R,

 First you need to clarify what you'd like to happen when the ID in B
 is not unique. What do you want the resulting dataframe to look
 like?

 Some possible answers involve using different options for merge() or
 using unique() to remove duplicates from B before merging. But
 at least to me, merge or retrieve the common ones isn't clear
 enough to be able to say which.

 Sarah

 On Wed, Aug 25, 2010 at 5:35 AM, Mangalani Peter Makananisa
 pmakanan...@sars.gov.za wrote:
 Dear R Gurus,



 I am currently working on the two dataset ( A and B), they both have the
 same fields:    ID , REGION, OFFICE, CSTART, CEND, NCYCLE, STATUS and
 CB.

 I want to merge the two data set by ID. The problem I have is that the
 in data A, the ID's are unique. However in the data set B, the ID's are
 not unique, thus some repeat themselves.



 How do I the merge or retrieve the common ones?

 Please advise.



 Kind Regards



 Peter



-- 
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] using an objects contents in a text string

2010-08-25 Thread Sarah Goslee
I'm a bit confused by your question, but you might just want paste:
paste(input.variable.name, data.txt, sep=_)

Sarah

On Wed, Aug 25, 2010 at 8:06 AM, josquint josqu...@unimelb.edu.au wrote:

 ALSO    I have had a play with cat() but have also not got this to
 work


 e.g.

 write(x, file = cat(input.variable.name  , file = , sep = _data.txt, ),
 ...)

 but this does not seem to work and I'm sure it is not the correct use of
 cat()

 Thanks.
 --



-- 
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] multi day R run crashed - is there a log

2010-08-25 Thread Martin Tomko

Hi Sarah,
thank you very much for your answer. I have been spitting things out on 
screen, but unfortunately I have not run it as a batch log, but in the 
interactive window, so when the GUI crashed I was left without trace... 
I guess I should google how to run it from a batch.


I should also explore using RData, I have ony been using csv files and 
tables so far, it seems that can also bring some added performance.


As the main output of my process is a matrix, I would really need to 
append to a matrix after each iteration. I have identified the 
write.table append parameter-based solution, but that would only append 
rows. Is there a way to slowly grow a matrix in both directions, 
meaning append columns as well (it is a big distance matrix).

Thanks a lot for your help,
Cheers
Martin

On 8/25/2010 1:56 PM, Sarah Goslee wrote:

Oh dear. No there isn't. I do a lot of very long runs, and have learned to write
out intermediate steps. Depending on what you are doing, saving the RData
file periodically may be appropriate, or writing out a csv file of
results so far,
or even just printing something to the output file (if you are using batch).

There may well be something more elegant, but these solutions have all worked
well for me in various situations.

Sarah

On Wed, Aug 25, 2010 at 3:46 AM, Martin Tomkomartin.to...@geo.uzh.ch  wrote:
   

Dear all,
I am using an R 2.10 installation on a Windows 203 server that I have no
control over. After a multi-day run I found that it was terminated/crashed.
Is there any log kept by R where I could see whether something/what
happened? The same process has been run beofre on a smaller dataset (also at
least a day of computing) without problems.

Thanks
Martin
 
   



--
Martin Tomko
Postdoctoral Research Assistant

Geographic Information Systems Division
Department of Geography
University of Zurich - Irchel
Winterthurerstr. 190
CH-8057 Zurich, Switzerland

email:  martin.to...@geo.uzh.ch
site:   http://www.geo.uzh.ch/~mtomko
mob:+41-788 629 558
tel:+41-44-6355256
fax:+41-44-6356848

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] multi day R run crashed - is there a log

2010-08-25 Thread Sarah Goslee
Hi,

On Wed, Aug 25, 2010 at 8:11 AM, Martin Tomko martin.to...@geo.uzh.ch wrote:
 Hi Sarah,
 thank you very much for your answer. I have been spitting things out on
 screen, but unfortunately I have not run it as a batch log, but in the
 interactive window, so when the GUI crashed I was left without trace... I
 guess I should google how to run it from a batch.

I have  no idea how to do that in Windows, but I'm sure it's possible. :)
That way the things written to the screen will be saved in a file instead.

 I should also explore using RData, I have ony been using csv files and
 tables so far, it seems that can also bring some added performance.

If you need to save *everything*, RData is what you get when you use save(),
or when you close a session and choose y to saving the data. It can be
read in using load(). That's one way to be able to pick up where you left
off.

 As the main output of my process is a matrix, I would really need to append
 to a matrix after each iteration. I have identified the write.table append
 parameter-based solution, but that would only append rows. Is there a way to
 slowly grow a matrix in both directions, meaning append columns as well
 (it is a big distance matrix).

AFAIK, you can't append columns that way because of the way text files are
written to disk. You'd need to rewrite the whole thing, or possibly write it out
in lower triangular format with NA values as padding (assuming it's a symmetric
distance).

Or for that matter, you could just write it out as a really long
vector, and turn
it back into a matrix later if you need to read the saved file in after a crash.

I'd recommend saving whatever variables are needed so that you can pick up
exactly where you left off, if possible. Much nicer to pick up 12 hours in than
to start over from the beginning.

Not R, but I just finished a 5-week batch job. You can bet that I put a lot of
thought into incremental save points and how to resume after an unexpected
halt!

Sarah

-- 
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] Estimate average standard deviation of mean of two dependent groups

2010-08-25 Thread Jokel Meyer
Dear R-experts!

I am currently running a meta-analysis with the help of the great metafor
package. However I have some difficulties setting up my raw data to enter it
to the meta-analysis models.
I have a group of subjects that have been measured in two continuous
variables (A  B). I have the information about the mean of the two
variables, the group size (n) and the standard deviations of the two
variables.
Now I would like to average both variables (A  B) and get the mean and
standard deviation of the merged variable (C).
As for the mean this would be quiet easy: I would just take the mean of mean
A and mean B to get the mean of C.
However for the standard deviation this seems more tricky as it is to assume
that standard deviations in A  B correlate. I assume (based on further
analysis) a correlation of r =0.5.
I found the formula to get the standard deviation of the SUM (not the mean)
of two variables:
SD=SQRT(SD_A^2 + SD_B^2 + 2*r*SD_A*SD_B)

with SD_B and SD_B being the standard deviation of A and B. And r*SD_A*SD_B
being the covariance of A and B.

Would this formula also be valid if I want to average (and not sum) my two
variables?

Many thanks for any help  best wishes,
Jokel

[[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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread Gabor Grothendieck
On Wed, Aug 25, 2010 at 7:43 AM, skan juanp...@gmail.com wrote:

  down vote  favorite


 Hello

 I have a zoo series. It lasts 10 years and its frequency is 15min.

 I'd like to get a new zoo series (or vector) with the same number of
 elements, whith each element equal to the first element of the day. That's,
 The first element everyday is repeated throughout the wole day.

 This is not same as aggregate(originalseries,as.Date,head,1) because this
 gives a vector with just one element for each day.

 cheers


Here are a few more solutions too:

library(zoo)
library(chron)
z - zoo(1:10, chron(0:9/5))

# aggregate / na.locf
z.ag - aggregate(z, as.Date, head, 1)
na.locf(z.ag, xout = time(z))

# duplicated / na.locf
z.na - ifelse.zoo(!duplicated(as.Date(time(z))), z, NA)
na.locf(z.na)

# ave - as before
zz - z
zz[] - ave(coredata(z), as.Date(time(z)), FUN = function(x) head(x, 1))
zz

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

2010-08-25 Thread Sarah Goslee
Almost. You'll need to handle the duplicate ID yourself since R has
no way of knowing which one(s) to change to NA. As I already suggested, you
can use unique() in conjunction with whatever logical rules you require
for choosing those values.

As I also already suggested, all.y and all.x are the options to merge()
that you need to consider.

 A - data.frame(ID = c(1,2,3), X = c('a','b','c'))
 B - data.frame(ID = c(1,2,2), Y = c('x','y','z'))
 merge(A, B, all.x=FALSE, all.y=TRUE)
  ID X Y
1  1 a x
2  2 b y
3  2 b z

Just think how much easier this process would have been if you had
provided a clear question with toy data and examples of what you'd
tried in your first question.

Sarah

On Wed, Aug 25, 2010 at 8:24 AM, Mangalani Peter Makananisa
pmakanan...@sars.gov.za wrote:
 A:
 ID X
 1  a
 2  b
 3  c


 B:
 ID Y
 1  x
 2  y
 2  z

 I would like to see something like this:

 Common = Merge(A,B)
 Common

 ID  X    Y
 1   a    x
 2   b    y
 2  N/A   z

 If it is possible,




-- 
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] modify a nomogram

2010-08-25 Thread Frank Harrell


Update to the rms package which is the version now being actively 
supported.  New features will not be added to Design.  The nomogram 
function in rms separates the plotting into a plot method for easier 
understanding.  You can control all axes - some experimentation can 
tell you if you can do what you want as I haven't tried exactly what 
you are doing.


Frank

Frank E Harrell Jr   Professor and ChairmanSchool of Medicine
 Department of Biostatistics   Vanderbilt University

On Wed, 25 Aug 2010, david dav wrote:


Hi,
I would like to emphasize (zoom) the zone of a nomogram where the
probability are  0.01
(nomogram built with nomogram, Design).
As a consequence, I don't need to draw the part of the Total points
axis with score  60 equivalent in my case to a linear predictor  4.5

- As far as I know, this is not possible with the arguments of the function.
- Changing the code of the function is beyond my abilities
-- can not even create a nomo.f function with the same body:
 body(nomo) - expression({
 conf.lp - match.arg(conf.lp) . rest of the function body
   this does not even work
-- I am not sure to find the piece of code responsible for
defining the axis

nomogram(logistic.fit, fun = plogis, fun.at = c(seq(0.01,1,by=.2)),
lmgp=.2, lp = T,maxscale = 100,   total.sep.page = T )

Thanks
David

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

2010-08-25 Thread Hadley Wickham
On Wed, Aug 25, 2010 at 6:05 AM, abotaha yaseen0...@gmail.com wrote:

 Woow, it is amazing,
 thank you very much.
 yes i forget to attach the dates, however, the dates in my case is every 16
 days.
 so how i can use 16 day interval instead of month in by option.

Here's one way using the lubridate package:

library(lubridate)
today() + days(16) * seq_len(50)
ymd(2008-01-01) + days(16) * seq_len(50)

-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

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


Re: [R] multi day R run crashed - is there a log

2010-08-25 Thread Martin Tomko

Hi ,
thanks, the lower.tri idea is I guess the best way. Will try that.
Cheers
Martin

On 8/25/2010 2:21 PM, Sarah Goslee wrote:

Hi,

On Wed, Aug 25, 2010 at 8:11 AM, Martin Tomkomartin.to...@geo.uzh.ch  wrote:
   

Hi Sarah,
thank you very much for your answer. I have been spitting things out on
screen, but unfortunately I have not run it as a batch log, but in the
interactive window, so when the GUI crashed I was left without trace... I
guess I should google how to run it from a batch.
 

I have  no idea how to do that in Windows, but I'm sure it's possible. :)
That way the things written to the screen will be saved in a file instead.

   

I should also explore using RData, I have ony been using csv files and
tables so far, it seems that can also bring some added performance.
 

If you need to save *everything*, RData is what you get when you use save(),
or when you close a session and choose y to saving the data. It can be
read in using load(). That's one way to be able to pick up where you left
off.

   

As the main output of my process is a matrix, I would really need to append
to a matrix after each iteration. I have identified the write.table append
parameter-based solution, but that would only append rows. Is there a way to
slowly grow a matrix in both directions, meaning append columns as well
(it is a big distance matrix).
 

AFAIK, you can't append columns that way because of the way text files are
written to disk. You'd need to rewrite the whole thing, or possibly write it out
in lower triangular format with NA values as padding (assuming it's a symmetric
distance).

Or for that matter, you could just write it out as a really long
vector, and turn
it back into a matrix later if you need to read the saved file in after a crash.

I'd recommend saving whatever variables are needed so that you can pick up
exactly where you left off, if possible. Much nicer to pick up 12 hours in than
to start over from the beginning.

Not R, but I just finished a 5-week batch job. You can bet that I put a lot of
thought into incremental save points and how to resume after an unexpected
halt!

Sarah

   



--
Martin Tomko
Postdoctoral Research Assistant

Geographic Information Systems Division
Department of Geography
University of Zurich - Irchel
Winterthurerstr. 190
CH-8057 Zurich, Switzerland

email:  martin.to...@geo.uzh.ch
site:   http://www.geo.uzh.ch/~mtomko
mob:+41-788 629 558
tel:+41-44-6355256
fax:+41-44-6356848

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


[R] multiple x factors using the sciplot package

2010-08-25 Thread Jordan Ouellette-Plante

Dear R community,

 

I am a beginner using the sciplot package to graph barplots. I would like to be 
able to graph two x factors (Sampling.year and Period). 

Sampling.year: 2006, 2007, 2008, 2009

Period: First, Second, Total

 

The parameter group is the different species I looked at. They can be seen in 
the legend.

The parameter response is the Average percentage cover category of these 
species.

 

I would like the graph so you see on the x axis the years (in a bigger font) 
and the period (in smaller font). It would look a bit like the barplot that can 
be seen at 
http://onertipaday.blogspot.com/2007/05/make-many-barplot-into-one-plot.html, 
but with the advantage of the sciplot package.

 

 

Thanks a lot for your help!

 

Jordan
  
[[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] percentage sign in expression

2010-08-25 Thread David Winsemius


On Aug 25, 2010, at 4:32 AM, e-letter wrote:


On 24/08/2010, David Winsemius dwinsem...@comcast.net wrote:


On Aug 24, 2010, at 9:37 AM, e-letter wrote:


Readers,

According to the documentation for the function 'plotmath' there  
is no

apparent possibility to add the percent sign (%) to a plot function,


Where did you see an assertion made???


Within R I entered the command:

?plotmath

Also accessed using:

help.start(browser=opera)

Navigating the web browser page:

packages
packages in /usr/lib/R/library
grdevices
plotmath

In the list headed 'syntax' and 'meaning' within the section  
'details'.



e.g.

plot(a[,1]~b[,2],ylab=expression(x~%),xlab=expression(z))





How to achieve this please?


Read the plotmath helo page more carefully. The section immediatedly
below the plotmath expressions points you to the use of the symbol()
expression-function and to the points help page where generation of
the available glyphs proceeds according to the advice on  
help(plotmath):



In my system the paragraph immediately after the list of features
(i.e. 'syntax','meaning') describes a note to TeX users. I cannot see
reference to 'symbol()'.


It's possible that my help page is different than yours. Right after  
the syntax/meaning description on mine (which is a Mac OSX system) is  
a paragraph:


The symbol font uses Adobe Symbol encoding so, for example, a lower  
case mu can be obtained either by the special symbol mu or by  
symbol(m). This provides access to symbols that have no special  
symbol name, for example, the universal, or forall, symbol is  
symbol(\042). To see what symbols are available in this way  
useTestChars(font=5) as given in the examples for points: some are  
only available on some devices.


(In this case I would be surprised if the help pages were different  
because this makes a cross-reference to the examples in points.  I am  
not surprised about cross-platform differences in descriptions of  
graphical devices  and would have included a caveat if I were  
corresponding on rhelp about such. I suppose the font issues could be  
platform specific so if you want to correct me on this point, I will  
try to file it away. I did, however, give you the code needed to to  
display Symbols and it sounds further on that it succeeded)





TestChars - function(sign=1, font=1, ...)

+ {
+if(font == 5) { sign - 1; r - c(32:126, 160:254)
+} else if (l10n_info()$MBCS) r - 32:126 else r - 32:255
+if (sign == -1) r - c(32:126, 160:255)
+par(pty=s)
+plot(c(-1,16), c(-1,16), type=n, xlab=, ylab=,
+ xaxs=i, yaxs=i)
+grid(17, 17, lty=1)
+for(i in r) try(points(i%%16, i%/%16, pch=sign*i,  
font=font,...))

+ }

TestChars(font=5)


Notice that the % sign is three characters to the right (i.e.
higher) of the forall symbol that is produced by the example code


I can't see 'forall' in the code above.


Gavin has already explained why you did not. The upside-down A (==  
universal or forall) was a useful reference point in the indexing,  
since it is only 3 glyphs away for the % symbol.





they offer. (The numbering proceeds from bottom up which confused me
at first.)


What numbering?


Actually I did not see any numbering either, which was why I remained  
confused about the location of the % symbol for several minutes.  
Perhaps I should have used the term indexing.





The documentation makes reference to the command:

demo(plotmath)

I applied this command and could not see an instruction to produce the
percent (%) symbol.


I don't think I suggested it would.

--

David Winsemius, MD
West Hartford, CT

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


Re: [R] Repeat the first day data through all the day. Zoo

2010-08-25 Thread skan

thanks
I'll try them, 

Why do you use the brackets in zz[]  ?
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Repeat-the-first-day-data-through-all-the-day-Zoo-tp2338069p2338266.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] how do R calculates the number of intervals between tick-marks

2010-08-25 Thread David Winsemius


On Aug 24, 2010, at 7:05 PM, Antonio Olinto wrote:


Hello,

I want to know how do R calculates the number of intervals between  
tick-marks in the y axis in a plot.


?axTicks # and then look at the other citations and the code as needed




I'm making a three y-axes plot and this information would help me a  
lot.


Thanks in advance.

Antonio Olinto


--

David Winsemius, MD
West Hartford, CT

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


[R] Comparing samples with widely different uncertainties

2010-08-25 Thread Sandy Small

Hi
This is probably more of a statistics question than a specific R
question, although I will be using R and need to know how to solve the
problem in R.

I have several sets of data (ejection fraction measurements) taken in
various ways from the same set of (~400) patients (so it is paired data).
For each individual measurement I can make an estimate of the percentage
uncertainty in the measurement.
Generally the measurements in data set A are higher but they have a
large uncertainty (~20%) while the measurements in data set Bare lower
but have a small uncertainty (~4%).
I believe, from the physiology, that the true value is likely to be
nearer the value of A than of B.
I need to show that, despite the uncertainties in the measurements
(which are not themselves normally distributed), there is (or is not) a
difference between the two groups, (a straight Wilcoxon signed ranks
test shows a difference but it cannot include that uncertainty data).

Can anybody suggest what I should be looking at? Is there a language
here that I don't know? How do I do it in R?
Many thanks for your help
Sandy

--

Sandy Small
Clinical Physicist
NHS Greater Glasgow and Clyde
and
NHS Forth Valley

Phone: 01412114592
E-mail: sandy.sm...@nhs.net




This message may contain confidential information. If yo...{{dropped:21}}

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread Gabor Grothendieck
On Wed, Aug 25, 2010 at 9:48 AM, skan juanp...@gmail.com wrote:

 thanks
 I'll try them,

 Why do you use the brackets in zz[]  ?

So it stays a zoo object with the same index.  We are only replacing
the data part.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 R calculates the number of intervals between tick-marks

2010-08-25 Thread Henrique Dallazuanna
Take a look at pretty function.

On Tue, Aug 24, 2010 at 8:05 PM, Antonio Olinto aolint...@bignet.com.brwrote:

 Hello,

 I want to know how do R calculates the number of intervals between
 tick-marks in the y axis in a plot.

 I'm making a three y-axes plot and this information would help me a lot.

 Thanks in advance.

 Antonio Olinto



 
 Webmail - iBCMG Internet
 http://www.ibcmg.com.br

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




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

[[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] Odp: Finding pairs

2010-08-25 Thread Petr PIKAL
Hi

well, I will add some explanation

r-help-boun...@r-project.org napsal dne 25.08.2010 11:24:38:

 Dear Mr Petr PIKAL
 After reading the R code provided by you, I realized that I would have 
never 
 figured out how this could have been done. I am going to re-read again 
and 
 again your code to understand the logic and the commands you have 
provided.
 Thanks again from the heart for your kind advice.
 Regards
 Mike
 
 --- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:
 
 From: Petr PIKAL petr.pi...@precheza.cz
 Subject: Re: [R] Odp:  Finding pairs
 To: Mike Rhodes mike_simpso...@yahoo.co.uk
 Cc: r-help@r-project.org
 Date: Wednesday, 25 August, 2010, 9:01
 
 Hm
 
 r-help-boun...@r-project.org napsal dne 25.08.2010 09:43:26:
 
  Dear Mr Petr Pikal
  
  I am extremely sorry for the manner I have raised the query. Actually 
 that was
  my first post to this R forum and in fact even I was also bit confused 

 while 
  drafting the query, for which I really owe sorry to all for consuming 
 the 
  precious time. Perhaps I will try to redraft my query in a better way 
as 
 follows. 
  
  I have two datasets A and B containing the names of branch offices 

 of a 
  particular bank say XYZ plc bank. The XYZ bank has number of main 
branch 
 
  offices (say Parent) and some small branch offices falling under the 
 purview 
  of these main branch offices (say Child).
  
  The datalist A and B consists of these main branch office names as 

 well as
  small branch office names. B is subset of A and these branch names are 

 coded. 
  Thus we have two datasets A and B as (again I am using only a
   portion of a large database just to have some idea)
  
  
  A B
  144  
what is here in B? Empty space?, 
  145   
  146   
  147  144
 
 How do you know that 144 from B relates to 147 in A? Is it according to 
 its positions? I.e. 4th item in B belongs to 4.th item in A?
 
  148  145  
  
  149  147
  151  148
  
  
  
  Now the branch 144 appears in A as well as in B and in B it is mapped 
 with 
  147. This means branch 147 comes under the purview of main branch 144. 

 Again 
  147 is controlling the branch 149 (since 147 also has appeared in B 
and 
 is 
  mapped with 149 of A).
  
  Similarly, branch 145 is controlling branch 148 which further controls 

  operations of bank branch 151 and like wise.
 
 Well as you did not say anything about structure of your data
 A-144:151
 B-144:148
 data.frame(A,B)
 A   B
 1 144  NA
 2 145  NA
 3 146  NA
 4 147 144
 5 148 145
 6 149 146
 7 150 147
 8 151 148
 DF-data.frame(A,B)

This was just making a data frame with 2 columns to have some data to play 
with

 main-DF$A[is.na(DF$B)]

Above are codes from A which are NA in B

 branch1-DF[!is.na(DF$B),]

Above is data frame of remaining codes (other than main)

 selected.branch1-branch1$A[branch1$B%in%main]

Above is codes from column A for which B column and main are the same

 branch2-branch1[!branch1$B%in%main,]

This is the rest of yet not selected rows

 selected.branch2-branch2$A[branch2$B%in%selected.branch1]

and this is selection of values from column A for which B column and 
selected.branch1 values are same.

But it works for this particular data, I am not sure how it behaves with 
duplicates and further issues. It also depends on how your data is 
organised.

And if you are in reading you could also go through setdiff, merge and 
maybe sqldf package and Rdata Import/export manual 

Regards
Petr


 
 and for cbinding your data which has uneven number of values see Jim 
 Holtman's answer to this
 
 How to cbind DF:s with differing number of rows?
 
 Regards
 Petr
 
 
  
  So in the end I need an output something like -
  
  Main Branch   Branch office1 Branch
   office2
  144 147  
   149
  145 148  
   151 

  146 NA
NA   
  
 
...
  
 
..
  
   
  I understand again I am not able to put forward my query properly. But 
I 
 must 
  thank all of you for giving a patient reading to my query and for 
 reverting 
  back earlier. Thanks once again.
  
  With warmest regards
  
  Mike 
  
  
  --- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:
  
  From: Petr PIKAL petr.pi...@precheza.cz
  Subject: Odp: [R] Finding
   pairs
  To: Mike Rhodes mike_simpso...@yahoo.co.uk
  Cc: r-help@r-project.org
  Date: Wednesday, 25 August, 2010, 6:39
  
  Hi
  
  without other details it is probably impossible to give you any 
 reasonable 
  advice. Do you 

[R] approxfun-problems (yleft and yright ignored)

2010-08-25 Thread Samuel Wuest
Dear all,

I have run into a problem when running some code implemented in the
Bioconductor panp-package (applied to my own expression data), whereby gene
expression values of known true negative probesets (x) are interpolated onto
present/absent p-values (y) between 0 and 1 using the *approxfun -
function*{stats}; when I have used R version 2.8, everything had
worked fine,
however, after updating to R 2.11.1., I got unexpected output (explained
below).

Please correct me here, but as far as I understand, the yleft and yright
arguments set the extreme values of the interpolated y-values in case the
input x-values (on whose approxfun is applied) fall outside range(x). So if
I run approxfun with yleft=1 and yright=0 with y-values between 0 and 1,
then I should never get any values higher than 1. However, this is not the
case, as this code-example illustrates:

 ### define the x-values used to construct the approxfun, basically these
are 2000 expression values ranging from ~ 3 to 7:
 xNeg - NegExprs[, 1]
 xNeg - sort(xNeg, decreasing = TRUE)

 ### generate 2000 y-values between 0 and 1:
 yNeg - seq(0, 1, 1/(length(xNeg) - 1))
 ### define yleft and yright as well as the rule to clarify what should
happen if input x-values lie outside range(x):
 interp - approxfun(xNeg, yNeg, yleft = 1, yright = 0, rule=2)
Warning message:
In approxfun(xNeg, yNeg, yleft = 1, yright = 0, rule = 2) :
  collapsing to unique 'x' values
 ### apply the approxfun to expression data that range from ~2.9 to 13.9
and can therefore lie outside range(xNeg):
  PV - sapply(AllExprs[, 1], interp)
 range(PV)
[1]0.000 6208.932
 summary(PV)
 Min.   1st Qu.Median  Mean   3rd Qu.  Max.
0.000e+00 0.000e+00 2.774e-03 1.299e+00 3.164e-01 6.209e+03

So the resulting output PV object contains data ranging from 0 to 6208, the
latter of which lies outside yleft and is not anywhere close to extreme
y-values that were used to set up the interp-function. This seems wrong to
me, and from what I understand, yleft and yright are simply ignored?

I have attached a few histograms that visualize the data distributions of
the objects I xNeg, yNeg, AllExprs[,1] (== input x-values) and PV (the
output), so that it is easier to make sense of the data structures...

Does anyone have an explanation for this or can tell me how to fix the
problem?

Thanks a million for any help, best, Sam

 sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-apple-darwin9.8.0

locale:
[1] en_IE.UTF-8/en_IE.UTF-8/C/C/en_IE.UTF-8/en_IE.UTF-8

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

other attached packages:
[1] panp_1.18.0   affy_1.26.1   Biobase_2.8.0

loaded via a namespace (and not attached):
[1] affyio_1.16.0 preprocessCore_1.10.0


-- 
-
Samuel Wuest
Smurfit Institute of Genetics
Trinity College Dublin
Dublin 2, Ireland
Phone: +353-1-896 2444
Web: http://www.tcd.ie/Genetics/wellmer-2/index.html
Email: wue...@tcd.ie
--
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] Comparing samples with widely different uncertainties

2010-08-25 Thread peter dalgaard

On Aug 25, 2010, at 3:57 PM, Sandy Small wrote:

 Hi
 This is probably more of a statistics question than a specific R
 question, although I will be using R and need to know how to solve the
 problem in R.
 
 I have several sets of data (ejection fraction measurements) taken in
 various ways from the same set of (~400) patients (so it is paired data).
 For each individual measurement I can make an estimate of the percentage
 uncertainty in the measurement.
 Generally the measurements in data set A are higher but they have a
 large uncertainty (~20%) while the measurements in data set Bare lower
 but have a small uncertainty (~4%).
 I believe, from the physiology, that the true value is likely to be
 nearer the value of A than of B.
 I need to show that, despite the uncertainties in the measurements
 (which are not themselves normally distributed), there is (or is not) a
 difference between the two groups, (a straight Wilcoxon signed ranks
 test shows a difference but it cannot include that uncertainty data).
 
 Can anybody suggest what I should be looking at? Is there a language
 here that I don't know? How do I do it in R?
 Many thanks for your help
 Sandy

Hm, well...

I don't think the issue is entirely well-defined, but let me try and give you 
some pointers:

For bivariate normal data (X,Y), the situation is that even if V(X) != V(Y) you 
still end up looking at X-Y if the question is whether the means are the same. 
It's sort of the only thing you _can_ do...

For non-normal data, it is not clear what the null hypothesis really is. The 
signed-rank test assumes that X-Y has a symmetric distribution, which is 
dubious if X is not symmetric and its variation dominates that of Y. You could 
also do a sign test and see if the differences has a median of zero (bear in 
mind that the median of a difference is different from the difference of the 
medians, but it could actually suffice.)

I'd probably start off with a simple plot of Y vs X and look for fan-out 
effects indicating that the variance depends on the mean. If it does, perhaps 
take logs or square roots and see if it makes the variances appear more stable 
and perhaps improve on the normality. Then maybe just do a paired t-test. 

-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk  Priv: pda...@gmail.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] Comparing samples with widely different uncertainties

2010-08-25 Thread Marc Schwartz
On Aug 25, 2010, at 8:57 AM, Sandy Small wrote:

 Hi
 This is probably more of a statistics question than a specific R
 question, although I will be using R and need to know how to solve the
 problem in R.
 
 I have several sets of data (ejection fraction measurements) taken in
 various ways from the same set of (~400) patients (so it is paired data).
 For each individual measurement I can make an estimate of the percentage
 uncertainty in the measurement.
 Generally the measurements in data set A are higher but they have a
 large uncertainty (~20%) while the measurements in data set Bare lower
 but have a small uncertainty (~4%).
 I believe, from the physiology, that the true value is likely to be
 nearer the value of A than of B.
 I need to show that, despite the uncertainties in the measurements
 (which are not themselves normally distributed), there is (or is not) a
 difference between the two groups, (a straight Wilcoxon signed ranks
 test shows a difference but it cannot include that uncertainty data).
 
 Can anybody suggest what I should be looking at? Is there a language
 here that I don't know? How do I do it in R?
 Many thanks for your help
 Sandy


The first place that I would start is at Martin Bland's page pertaining to the 
design and analysis of measurement studies:

  http://www-users.york.ac.uk/~mb55/meas/meas.htm

The papers he co-authored with Doug Altman are the go to resource for this 
domain.

HTH,

Marc Schwartz

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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] SEM : Warning : Could not compute QR decomposition of Hessian

2010-08-25 Thread Anne Mimet

Hi useRs,

I'm trying for the first time to use a sem. The model finally runs,  
but gives a warning saying :
In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names  
= vars,  :   Could not compute QR decomposition of Hessian.

Optimization probably did not converge. 

I found in R-help some posts on this warning, but my attemps to modify  
the code didn't change the warning message (i tried to give an error  
of 1 to the latente variables). I can't figure what the problem is.

Here is the code :


tab-read.table(F:/Mes documents/stats/sem/donnees_corr.txt,  
header=T, sep=,na.strings = NA)


tab[,46]-as.factor(tab[,46])
tab[,24]-as.factor(tab[,24])
tab[,40]-as.factor(tab[,40])

fct_cor-hetcor(tab, ML=T)
cor_tab- fct_cor$correlations

moment_tab-read.moments(diag=F, names=c('c1','c2', 'c3','c4','c5',  
'c6','c7', 'c8', 'c9',  'ind_plando', 'long_sup15', 'long_inf15',  
'pente', 'est', 'sud','ouest', 'nord' ,'reg_hydriq', 'prof_sol',  
'pierro', 'efferv', 'struct','drainage','texture', 'route1_pond',  
'route2_pond',
'pourcactif', 'tx_chomage', 'pourcagric', 'pourc_jeunes', 'pop99',  
'rev_imp_foyer','eq_CONC', 'eq_sante', 'eq_edu', 'sold_nat',
'sold_mig', 'tps_dom_emp','TXEMPLOI','ORIECO','dist_paris','axe1',  
'axe2', 'axe3', 'dist_protect','urbanisation','pays_incli','pays_alti'))

# after comes the moment matrix (triangular)

ram_tab-specify.model()
type_paysage-pays_alti,NA,1
type_paysage-pays_incli, pays2, NA
pedo-reg_hydriq, NA, 1
pedo-prof_sol, ped8, NA
pedo-pierro, ped9, NA
pedo-efferv, ped10, NA
pedo-struct, ped11, NA
pedo-drainage, ped12, NA
pedo-texture, ped13, NA
adj_99-c1, NA,1
adj_99-c2, adj2,NA
adj_99-c3, adj3,NA
adj_99-c4, adj4,NA
adj_99-c5, adj5,NA
adj_99-c6, adj6,NA
adj_99-c7, adj7,NA
adj_99-c8, adj8,NA
adj_99-c9, adj9,NA
etat_hexa-axe1, NA, 1
etat_hexa-axe2, et2, NA
etat_hexa-axe3, et3, NA
socioBV-sold_mig, BV1, NA
socioBV-sold_nat, BV2, NA
socioBV-TXEMPLOI, BV3, NA
socioBV-ORIECO, BV4, NA
socioBV-tps_dom_emp, NA, 1
eqBV-eq_CONC, NA, 1
eqBV-eq_sante, eq2, NA
eqBV-eq_edu, eq3, NA
socio_com-pourcactif , NA, 1
socio_com-tx_chomage, com2, NA
socio_com-pourcagric, com3, NA
socio_com-pourc_jeunes, com4, NA
socio_com-pop99, com5, NA
socio_com-rev_imp_foyer, com7, NA
access_hexa-route1_pond, NA, 1
access_hexa-route2_pond, acc2, NA
hydro-ind_plando, NA, 1
hydro-long_sup15, eau2, NA
hydro-long_inf15, eau3, NA
topog-pente, NA, 1
topog-est, top2, NA
topog-sud, top3, NA
topog-nord, top4, NA
topog-ouest, top5, NA
dist_protect- urbanisation, cor1,NA
dist_protect- adj_99, cor2, NA
dist_protect- etat_hexa, cor3, NA
topog- urbanisation, cor4, NA
topog- adj_99, cor5, NA
topog- etat_hexa, cor6, NA
topog- access_hexa, cor7, NA
topog-hydro, cor8, NA
topog-pedo, cor9, NA
pedo- urbanisation, cor10, NA
pedo- adj_99, cor11, NA
pedo- etat_hexa, cor12, NA
pedo-hydro, cor1, NA
hydro- urbanisation, cor13, NA
hydro- adj_99, cor14, NA
hydro- etat_hexa, cor15, NA
access_hexa- urbanisation, cor16, NA
access_hexa- etat_hexa, cor17, NA
socio_com- etat_hexa, cor18, NA
socio_com- adj_99, cor19, NA
socio_com- urbanisation, cor20, NA
dist_paris- socio_com, cor21, NA
dist_paris- access_hexa, cor22, NA
dist_paris- adj_99, cor23, NA
dist_paris- etat_hexa, cor24, NA
dist_paris- urbanisation, cor25, NA
dist_paris- socioBV, cor26, NA
socioBV- eqBV, cor27, NA
socioBV- urbanisation, cor28, NA
socioBV- adj_99, cor29, NA
socioBV- etat_hexa, cor30, NA
eqBV- etat_hexa, cor31, NA
eqBV- adj_99, cor32, NA
eqBV- urbanisation, cor33, NA
etat_hexa- urbanisation, cor34, NA
etat_hexa- adj_99, cor35, NA
adj_99- urbanisation, cor36, NA
type_paysage- urbanisation, cor37, NA
type_paysage- adj_99, cor38, NA
type_paysage- etat_hexa, cor39, NA
dist_paris-dist_paris, auto1, NA
dist_protect-dist_protect, auto2, NA
c1 - c1, auto4, NA
c2 - c2 , auto5, NA
c3 -  c3 , auto6, NA
c4 -  c4 , auto7, NA
c5 -  c5 , auto8, NA
c6 -  c6 , auto9, NA
c7 -  c7 , auto10, NA
c8 - c8 , auto11, NA
c9 - c9 , auto12, NA
ind_plando -  ind_plando, auto13, NA
long_sup15 - long_sup15 , auto14, NA
long_inf15 -  long_inf15  , auto15, NA
pente - pente , auto16, NA
est-  est , auto17, NA
sud -  sud , auto18, NA
ouest- ouest , auto19, NA
nord -  nord , auto20, NA
reg_hydriq - reg_hydriq , auto21, NA
prof_sol- prof_sol , auto22, NA
pierro - pierro , auto23, NA
efferv -  efferv , auto24, NA
struct - struct , auto25, NA
drainage - drainage, auto26, NA
texture -  texture , auto27, NA
route1_pond  -route1_pond  , auto30, NA
route2_pond - route2_pond , auto31, NA
pourcactif  -  pourcactif , auto32, NA
tx_chomage -  tx_chomage, auto33, NA
pourcagric - pourcagric , auto34, NA
pourc_jeunes- pourc_jeunes  , auto36, NA
pop99- pop99  , auto36, NA
rev_imp_foyer -   rev_imp_foyer , auto38, NA
eq_CONC - eq_CONC, auto39, NA
eq_sante -eq_sante , auto40, NA
eq_edu-eq_edu   , auto41, NA
sold_nat-  sold_nat , auto42, NA
sold_mig  -  sold_mig  , auto43, NA
tps_dom_emp  -  tps_dom_emp , auto44, NA
TXEMPLOI  - TXEMPLOI  , auto45, NA
ORIECO  - ORIECO  , auto46, NA
axe1  -  axe1 , auto47, NA

[R] frequency, count rows, data for heat map

2010-08-25 Thread rtsweeney

Hi all, 
I have read posts of heat map creation but I am one step prior --
Here is what I am trying to do and wonder if you have any tips?
We are trying to map sequence reads from tumors to viral genomes.

Example input file :
111 abc
111 sdf
111 xyz
1079   abc
1079   xyz
1079   xyz
5576   abc
5576   sdf
5576   sdf

How may xyz's are there for 1079 and 111? How many abc's, etc?
How many times did reads from sample (1079) align to virus xyz. 
In some cases there are thousands per virus in a give sample, sometimes one.
The original file (two columns by tens of thousands of rows; 20 MB) is
text file (tab delimited).

Output file:
 abc  sdf  xyz
111 1  1 1
1079   1  0 2
5576   1  2 0

Or, other ways to generate this data so I can then use it for heat map
creation? 

Thanks for any help you may have, 

rtsweeney
palo alto, ca
-- 
View this message in context: 
http://r.789695.n4.nabble.com/frequency-count-rows-data-for-heat-map-tp2338363p2338363.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] lattice help required

2010-08-25 Thread Kay Cichini

hello,

i want to stack two lattice plots beneath each other using one x-axis and
sharing the same text-panels,
like:
#
library(lattice)

y1 - rnorm(100,100,10)
y2 - rnorm(100,10,1)
facs-expand.grid(Sites=rep(c(Site I,Site II),25),Treatment=c(A,B))
pl1-dotplot(y1 ~ facs$Treatment|facs$Sites,
 scales = list(x = list(rot = 90, tck=c(1,0
pl2-dotplot(y2 ~ facs$Treatment|facs$Sites,
 scales = list(x = list(rot = 90, tck=c(1,0

print(pl1, split=c(1,2,1,2), more=TRUE)
print(pl2, split=c(1,1,1,2))
#

but as said, ideally the plots should be stacked with only the lower plot
giving the x-axis annotation 
and only the upper plot with text-panels.

thanks a lot,
kay

-

Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck


-- 
View this message in context: 
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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] lmer() causes segfault

2010-08-25 Thread Bertolt Meyer

Dennis,

just wow. Thank you so much. I knew it was something trivial - in this  
case the variable type of the of the grouping variables. However,  
something as trivial as this should not throw a segfault IMHO. I tried  
subscribing to R-sig-mixed this morning, but the corresponding mail  
server at the ETH's stats department seems to be down. And thank you  
so much for changing the model, that is a great new starting point.  
Can you recommend a good book that deals with multilevel models in  
lmer() that include longitudinal data? I was not aware of the  
difference between scalar random effects and random slopes and would  
like to read up on that.


Again, thanks a lot.
Regards,
Bertolt


Am 25.08.2010 um 13:47 schrieb Dennis Murphy:


Hi:

Let's start with the data:

 str(test.data)
'data.frame':   100 obs. of  4 variables:
 $ StudentID: num  17370 17370 17370 17370 17379 ...
 $ GroupID  : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Time : num  1 2 3 4 1 2 3 4 1 2 ...
 $ Score: num  76.8 81.8 89.8 92.8 75.9 ...

Both StudentID and GroupID are numeric; in the model, they would be  
treated as continuous covariates rather than factors, so we need to  
convert:


test.data$StudentID - factor(test.data$StudentID)
test.data$GroupID - factor(test.data$GroupID)

Secondly, I believe there are some flaws in your model. After  
converting your variables to factors, I ran


library(lme4)
mlmoded1.lmer - lmer(Score ~ Time + (Time | GroupID/StudentID),  
data = test.data)


You have two groups, so they should be treated as a fixed effect -  
more specifically, as a fixed blocking factor. The StudentIDs are  
certainly nested within GroupID, and Time is measured on each  
StudentID, so it is a repeated measures factor. The output of this  
model is

 mlmoded1.lmer
Linear mixed model fit by REML
Formula: Score ~ Time + (Time | GroupID/StudentID)
   Data: test.data
   AIC   BIC logLik deviance REMLdev
 393.1 416.5 -187.5376.9   375.1
Random effects:
 GroupsNameVariance  Std.Dev. Corr
 StudentID:GroupID (Intercept)  0.504131 0.71002
   Time 0.083406 0.28880  1.000
 GroupID   (Intercept) 12.809567 3.57905
   Time 3.897041 1.97409  -1.000
 Residual   1.444532 1.20189
Number of obs: 100, groups: StudentID:GroupID, 25; GroupID, 2

Fixed effects:
Estimate Std. Error t value
(Intercept)   72.803  2.552  28.530
Time   4.474  1.401   3.193

Correlation of Fixed Effects:
 (Intr)
Time -0.994

The high correlations among the random effects and then among the  
fixed effects suggests that the model specification may be a bit off.


The above model fits random slopes to GroupIDs and StudentIDs, along  
with random intercepts, but GroupID is a between-subject effect and  
should be at the top level. Time is a within-subject effect and  
StudentIDs are the observational units. I modified the model to  
provide fixed effects for GroupIDs, scalar random effects for  
StudentIDs and random slopes for StudentIDs.


 mod3 - lmer(Score ~ 1 + GroupID + Time + (1 | StudentID) +
+  (0 + Time | StudentID), data = test.data)
 mod3
Linear mixed model fit by REML
Formula: Score ~ 1 + GroupID + Time + (1 | StudentID) + (0 + Time |  
StudentID)

   Data: test.data
   AIC   BIC logLik deviance REMLdev
 430.9 446.5 -209.4418.4   418.9
Random effects:
 GroupsNameVariance   Std.Dev.
 StudentID (Intercept) 4.2186e-13 6.4951e-07
 StudentID Time1.8380e+00 1.3557e+00
 Residual  1.6301e+00 1.2768e+00
Number of obs: 100, groups: StudentID, 25

Fixed effects:
Estimate Std. Error t value
(Intercept)  70.7705 0.4204  168.33
GroupID2  4.0248 0.58546.88
Time  4.5292 0.2942   15.39

Correlation of Fixed Effects:
 (Intr) GrpID2
GroupID2 -0.668
Time -0.264  0.000

I didn't check the quality of the fit, but on the surface it seems  
to be more stable, FWIW. Perhaps one could also add a term (GroupID  
| StudentID), but I don't know offhand if that would make any sense.  
Another issue to consider is whether to fit by REML or ML, but that  
is secondary to getting the form of the model equation right. I  
don't claim this as a final model, but rather a 're-starting point'.  
It may well be in need of improvement, so comments are welcome.


The confusion between subjects nested in time or vice versa has  
occurred several times this week with respect to repeated measures/ 
longitudinal models using lmer(), so perhaps it merits a comment:  
subjects/experimental units are NOT nested in time. Measurements  
taken on an individual at several time points *entails* that time be  
nested within subject. Just saying...


This discussion may be better continued on the R-sig-mixed list, so  
I've cc-ed to that group as well.


HTH,
Dennis

On Wed, Aug 25, 2010 at 1:27 AM, Bertolt Meyer  
bme...@sozpsy.uzh.ch wrote:

Ben Bolker bbolker at gmail.com 

Re: [R] frequency, count rows, data for heat map

2010-08-25 Thread Jan van der Laan
Your problem is not completely clear to me, but perhaps something like

data - data.frame(
   a = rep(c(1,2), each=10),
   b = rep(c('a', 'b', 'c', 'd'), 5))
library(plyr)
daply(data, a ~ b, nrow)

does what you need.

Regards,
Jan

On Wed, Aug 25, 2010 at 4:53 PM, rtsweeney tripswee...@gmail.com wrote:

 Hi all,
 I have read posts of heat map creation but I am one step prior --
 Here is what I am trying to do and wonder if you have any tips?
 We are trying to map sequence reads from tumors to viral genomes.

 Example input file :
 111     abc
 111     sdf
 111     xyz
 1079   abc
 1079   xyz
 1079   xyz
 5576   abc
 5576   sdf
 5576   sdf

 How may xyz's are there for 1079 and 111? How many abc's, etc?
 How many times did reads from sample (1079) align to virus xyz.
 In some cases there are thousands per virus in a give sample, sometimes one.
 The original file (two columns by tens of thousands of rows; 20 MB) is
 text file (tab delimited).

 Output file:
         abc  sdf  xyz
 111     1      1     1
 1079   1      0     2
 5576   1      2     0

 Or, other ways to generate this data so I can then use it for heat map
 creation?

 Thanks for any help you may have,

 rtsweeney
 palo alto, ca
 --
 View this message in context: 
 http://r.789695.n4.nabble.com/frequency-count-rows-data-for-heat-map-tp2338363p2338363.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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread skan

# duplicated / na.locf   doesn't work
it says Error in fix.by(by.x, x) : 'by' must specify valid column(s)

if I use ifelse instead of ifelse.zoo it works but it gives me a non zoo
vector.
Myabe is because my zoo version is older.

cheers
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Repeat-the-first-day-data-through-all-the-day-Zoo-tp2338069p2338409.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] Estimate average standard deviation of mean of two dependent groups

2010-08-25 Thread Joshua Wiley
Hi Jokel,

If I remember correctly

1) The variance of the sum of two variables X and Y is:
   var(X) + var(Y) + (2 * cov(X, Y))

   If you create some sample data, you can verify this by:
   var((X + Y))

2) The variance of a random variable (X) multiplied by some constant (K)
is equal to the the variance of X multiplied by the square of K.  So:

var(X * K) = var(X) * K^2

I have never seen these combined, but I imagine they could be.
Let Z = X + Y
The mean of X and Y is (X + Y)/2 = Z/2
From the formula above (1), the variance of Z is:

var(X) + var(Y) + (2 * cov(X, Y))

Because Z/2 = Z * 0.5, the variance of Z * 0.5 is given by:

var(Z) * (0.5^2)

and substituting, the variance of the mean of X and Y is:

(var(X) + var(Y) + (2 * cov(X, Y))) * (0.5^2)

This held up under several empirical tests I did where I had actual
data for X and Y.

***Discalimer
I am piecing this together, and I am far from an expert on the
subject.  I do not know if it is actually true, and there may be
additional assumptions.

I tested this using:

myfun - function() {
  x - rnorm(100)
  y - rnorm(100)
  var.x - var(x)
  var.y - var(y)
  cov.xy - cov(x, y)
  calc.var.xplusy - var.x + var.y + 2*cov.xy
  var.meanxy - calc.var.xplusy * (.5^2)
  empirical.var.meanxy - var( (x + y)/2 )
  output - all.equal(var.meanxy, empirical.var.meanxy)
  return(output)
}

temp - vector(logical, 1000)
for(i in 1:1000) {temp[i] - myfun()}
all(temp)

HTH,

Josh

On Wed, Aug 25, 2010 at 5:29 AM, Jokel Meyer jokel.me...@googlemail.com wrote:
 Dear R-experts!

 I am currently running a meta-analysis with the help of the great metafor
 package. However I have some difficulties setting up my raw data to enter it
 to the meta-analysis models.
 I have a group of subjects that have been measured in two continuous
 variables (A  B). I have the information about the mean of the two
 variables, the group size (n) and the standard deviations of the two
 variables.
 Now I would like to average both variables (A  B) and get the mean and
 standard deviation of the merged variable (C).
 As for the mean this would be quiet easy: I would just take the mean of mean
 A and mean B to get the mean of C.
 However for the standard deviation this seems more tricky as it is to assume
 that standard deviations in A  B correlate. I assume (based on further
 analysis) a correlation of r =0.5.
 I found the formula to get the standard deviation of the SUM (not the mean)
 of two variables:
 SD=SQRT(SD_A^2 + SD_B^2 + 2*r*SD_A*SD_B)

 with SD_B and SD_B being the standard deviation of A and B. And r*SD_A*SD_B
 being the covariance of A and B.

 Would this formula also be valid if I want to average (and not sum) my two
 variables?

 Many thanks for any help  best wishes,
 Jokel

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




-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.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] lattice help required

2010-08-25 Thread RICHARD M. HEIBERGER
Kay,

doe this do what you want?


dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(x = list(rot = 90, tck=c(1,0))),
ylab=c(y1, y2),
xlab=c(Site 1, Site 2),
strip=FALSE)


On Wed, Aug 25, 2010 at 11:04 AM, Kay Cichini kay.cich...@uibk.ac.atwrote:


 hello,

 i want to stack two lattice plots beneath each other using one x-axis and
 sharing the same text-panels,
 like:
 #
 library(lattice)

 y1 - rnorm(100,100,10)
 y2 - rnorm(100,10,1)
 facs-expand.grid(Sites=rep(c(Site I,Site II),25),Treatment=c(A,B))
 pl1-dotplot(y1 ~ facs$Treatment|facs$Sites,
 scales = list(x = list(rot = 90, tck=c(1,0
 pl2-dotplot(y2 ~ facs$Treatment|facs$Sites,
 scales = list(x = list(rot = 90, tck=c(1,0

 print(pl1, split=c(1,2,1,2), more=TRUE)
 print(pl2, split=c(1,1,1,2))
 #

 but as said, ideally the plots should be stacked with only the lower plot
 giving the x-axis annotation
 and only the upper plot with text-panels.

 thanks a lot,
 kay

 -
 
 Kay Cichini
 Postgraduate student
 Institute of Botany
 Univ. of Innsbruck
 

 --
 View this message in context:
 http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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.htmlhttp://www.r-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


[[alternative HTML version deleted]]

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


Re: [R] Draw a perpendicular line?

2010-08-25 Thread William Revelle

At 3:04 PM -0700 8/23/10, CZ wrote:

Hi,

I am trying to draw a perpendicular line from a point to two points.
Mathematically I know how to do it, but to program it, I encounter some
problem and hope can get help.  Thanks.

I have points, A, B and C.  I calculate the slope and intercept for line
drawn between A and B.
I am trying to check whether I can draw a perpendicular line from C to line
AB and get the x,y value for the point D at the intersection.

Assume I get the slope of the perpendicular line, I will have my point (D)
using variable x and y which is potentially on line AB.   My idea was using
|AC|*|AC| = |AD|*|AD|+ |CD|*|CD|.  I don't know what function I may need to
call to calculate the values for point D (uniroot?).



This is easier than you think.
Think of the x,y coordinates of each point :
Then, the slope is slope = rise/run =  (By- Ay)/(Bx- Ax)
The Dx coordinate = Cx and the Dy = (Dx - Ax) * slope
Then, to draw the line segment from C to D
lines(C,D)

In R:

A - c(2,4)
B - c(4,1)
C - c(8,10)
slope -( C[2]- A[2])/(C[1]-A[1])  #rise/run
D - c(B[1],(B[1]-A[1])*slope + A[2])  # find D

my.data - rbind(A,B,C,D)
colnames(my.data) - c(X,Y)
my.data#show it
plot(my.data,type=n)   #graph it without the points
text(my.data,rownames(my.data))  #put the points in
segments(A[1],A[2],C[1],C[2])   #draw the line segments
segments(B[1],B[2],D[1],D[2])

Bill




Thank you.



--
View this message in context: 
http://r.789695.n4.nabble.com/Draw-a-perpendicular-line-tp2335882p2335882.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] Change value of a slot of an S4 object within a method.

2010-08-25 Thread Joris Meys
Dear all,

I have an S4 class with a slot extra which is a list. I want to be
able to add an element called name to that list, containing the
object value (which can be a vector, a dataframe or another S4
object)

Obviously

setMethod(add.extra,signature=c(PM10Data,character,vector),
  function(object,name,value){
 obj...@extra[[name]] - value
  }
)

Contrary to what I would expect, the line :
eval(eval(substitute(expression(obj...@extra[[name]] - value

gives the error :
Error in obj...@extra[[name]] - value : object 'object' not found

Substitute apparently doesn't work any more in S4 methods...

 I found a work-around by calling the initializer every time, but this
influences the performance. Returning the object is also not an
option, as I'd have to remember to assign that one each time to the
original name.

Basically I'm trying to do some call by reference with S4, but don't
see how I should do that. How would I be able to tackle this problem
in an efficient and elegant way?

Thank you in advance
Cheers
Joris

-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
joris.m...@ugent.be
---
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php

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


Re: [R] multiple x factors using the sciplot package

2010-08-25 Thread Dennis Murphy
Hi:

You can probably do what you want in either ggplot2 or lattice, but I would
recommend at least a couple different approaches:

(1) Plot individual bar charts by combinations of year and period.
 This is easy to do in both ggplot2 and lattice: in ggplot2, one
 would use  geom_bar(x)  + facet_grid(year ~ period), where
 x is the count variable. The same thing in lattice would be
 barchart( ~ x | year * period, data = mydata). Suitable options
 can be added to either plot for better presentation. All the plots
 are on the same horizontal and vertical scales so that visual
 comparisons across the rows and columns are easy to make.
(2) Use a Cleveland dot chart. This is much less wasteful of ink
  than bar charts, and multiple groups can be compared without
  too much difficulty. See dotplot() in lattice for details. The
  VADeaths example is likely the best one to emulate for your
  problem. Both the input data structure and plot code matter.
(3) If your purpose is to compare the two groups with respect
  to a count variable, a third option is to consider a mosaic plot,
  which can be thought of as a 'visual chi-square test'. These
  can be found in package vcd, which contains a nice 40+
  page vignette to show you in detail how to create a suitable
  mosaic plot (and how to interpret it).

I prefer to see R used to promote good graphics and sound statistical
principles. In that spirit, performing a two-factor comparison with count
data is easier, and quite possibly more instructive, with a 2D visual
metaphor (which all three of the suggestions above have in common) than with
the 1D metaphor of multiple dodged/side-by-side bar charts. Of course, this
may be against the standard practice in your field, in which case you have a
decision to make.

HTH,
Dennis

On Wed, Aug 25, 2010 at 6:16 AM, Jordan Ouellette-Plante 
jordan_opla...@hotmail.com wrote:


 Dear R community,



 I am a beginner using the sciplot package to graph barplots. I would like
 to be able to graph two x factors (Sampling.year and Period).

 Sampling.year: 2006, 2007, 2008, 2009

 Period: First, Second, Total



 The parameter group is the different species I looked at. They can be
 seen in the legend.

 The parameter response is the Average percentage cover category of these
 species.



 I would like the graph so you see on the x axis the years (in a bigger
 font) and the period (in smaller font). It would look a bit like the barplot
 that can be seen at
 http://onertipaday.blogspot.com/2007/05/make-many-barplot-into-one-plot.html,
 but with the advantage of the sciplot package.





 Thanks a lot for your help!



 Jordan

[[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] Compiling Fortran for R : .so: mach-o, but wrong architecture

2010-08-25 Thread holstius

Hi Marie, this link may be helpful if you want to build it for both
architectures, 32-bit and 64-bit. I struggled with the same problem not too
long ago. 

 
http://cran.rakanu.com/bin/macosx/RMacOSX-FAQ.html#Building-universal-package  

I ended up writing a very simple Makefile to accompany [my version of] the
bar.f code. There is almost certainly a better way to do it. Maybe with
Makevars?

David 
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Compiling-Fortran-for-R-so-mach-o-but-wrong-architecture-tp2337198p2338507.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] Change value of a slot of an S4 object within a method.

2010-08-25 Thread Steve Lianoglou
Hi Joris,

On Wed, Aug 25, 2010 at 11:56 AM, Joris Meys jorism...@gmail.com wrote:
 Dear all,

 I have an S4 class with a slot extra which is a list. I want to be
 able to add an element called name to that list, containing the
 object value (which can be a vector, a dataframe or another S4
 object)

 Obviously

 setMethod(add.extra,signature=c(PM10Data,character,vector),
  function(object,name,value){
             obj...@extra[[name]] - value
  }
 )

 Contrary to what I would expect, the line :
 eval(eval(substitute(expression(obj...@extra[[name]] - value

 gives the error :
 Error in obj...@extra[[name]] - value : object 'object' not found

 Substitute apparently doesn't work any more in S4 methods...

  I found a work-around by calling the initializer every time, but this
 influences the performance. Returning the object is also not an
 option, as I'd have to remember to assign that one each time to the
 original name.

 Basically I'm trying to do some call by reference with S4, but don't
 see how I should do that. How would I be able to tackle this problem
 in an efficient and elegant way?

In lots of my own S4 classes I define a slot called .cache which is
an environment for this exact purpose.

Using this solution for your scenario might look something like this:

setMethod(add.extra,signature=c(PM10Data,character,vector),
function(object, name, value) {
  obj...@.cache$extra[[name]] - value
})

I'm not sure what your entire problem looks like, but to get your
extra list, or a value form it, you could:

setMethod(extra, signature=PM10Data,
function(object, name=NULL) {
  if (!is.null(name)) {
obj...@.cache$extra[[name]]
  } else {
obj...@.cache$extra
})

... or something like that.

The last thing you have to be careful of is that you nave to make sure
that each new(PM10Data) object you have initializes its *own* cache:

setClass(PM10Data, representation(..., .cache='environment'))
setMethod(initialize, PM10Data,
function(.Object, ..., .cache=new.env()) {
  callNextMethod(.Object, .cache=.cache, ...)
})

Does that make sense?

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

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


Re: [R] Odp: Finding pairs

2010-08-25 Thread Dennis Murphy
Hi:

I'm just ideating here (think IBM commercial...) but perhaps a graphical
model approach might be worth looking into. It seems to me that Mr. Rhodes
is looking for clusters of banks that are under the same ownership umbrella.
That information is not directly available in a single variable, but can
evidently be inferred from the matches between the two variables: B[i]
controls A[i] if B[i] is nonempty. In the bank 144 - 147 - 149 example,
149 controls 147 and 147 controls 144, so it appears that some transitive
relation holds among the set of matches as well. (Why is PacMan going
through my head? :)  I know next to nothing about graphical models, but I'm
thinking about igraph and some of the tools in the statnet bundle to tackle
this problem. Does that make sense to anyone? Alternatives?

FWIW,
Dennis

On Wed, Aug 25, 2010 at 2:24 AM, Mike Rhodes mike_simpso...@yahoo.co.ukwrote:

 Dear Mr Petr PIKAL
 After reading the R code provided by you, I realized that I would have
 never figured out how this could have been done. I am going to re-read again
 and again your code to understand the logic and the commands you have
 provided.
 Thanks again from the heart for your kind advice.
 Regards
 Mike

 --- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:

 From: Petr PIKAL petr.pi...@precheza.cz
 Subject: Re: [R] Odp:  Finding pairs
 To: Mike Rhodes mike_simpso...@yahoo.co.uk
 Cc: r-help@r-project.org
 Date: Wednesday, 25 August, 2010, 9:01

 Hm

 r-help-boun...@r-project.org napsal dne 25.08.2010 09:43:26:

  Dear Mr Petr Pikal
 
  I am extremely sorry for the manner I have raised the query. Actually
 that was
  my first post to this R forum and in fact even I was also bit confused
 while
  drafting the query, for which I really owe sorry to all for consuming
 the
  precious time. Perhaps I will try to redraft my query in a better way as
 follows.
 
  I have two datasets A and B containing the names of branch offices
 of a
  particular bank say XYZ plc bank. The XYZ bank has number of main branch

  offices (say Parent) and some small branch offices falling under the
 purview
  of these main branch offices (say Child).
 
  The datalist A and B consists of these main branch office names as
 well as
  small branch office names. B is subset of A and these branch names are
 coded.
  Thus we have two datasets A and B as (again I am using only a
   portion of a large database just to have some idea)
 
 
  A B
  144
what is here in B? Empty space?,
  145
  146
  147  144

 How do you know that 144 from B relates to 147 in A? Is it according to
 its positions? I.e. 4th item in B belongs to 4.th item in A?

  148  145
 
  149  147
  151  148
 
 
 
  Now the branch 144 appears in A as well as in B and in B it is mapped
 with
  147. This means branch 147 comes under the purview of main branch 144.
 Again
  147 is controlling the branch 149 (since 147 also has appeared in B and
 is
  mapped with 149 of A).
 
  Similarly, branch 145 is controlling branch 148 which further controls
  operations of bank branch 151 and like wise.

 Well as you did not say anything about structure of your data
 A-144:151
 B-144:148
 data.frame(A,B)
 A   B
 1 144  NA
 2 145  NA
 3 146  NA
 4 147 144
 5 148 145
 6 149 146
 7 150 147
 8 151 148
 DF-data.frame(A,B)
 main-DF$A[is.na(DF$B)]
 branch1-DF[!is.na(DF$B),]
 selected.branch1-branch1$A[branch1$B%in%main]
 branch2-branch1[!branch1$B%in%main,]
 selected.branch2-branch2$A[branch2$B%in%selected.branch1]

 and for cbinding your data which has uneven number of values see Jim
 Holtman's answer to this

 How to cbind DF:s with differing number of rows?

 Regards
 Petr


 
  So in the end I need an output something like -
 
  Main Branch   Branch office1 Branch
   office2
  144 147 149
  145 148 151

  146 NA
NA
 

 ...
 

 ..
 
 
  I understand again I am not able to put forward my query properly. But I
 must
  thank all of you for giving a patient reading to my query and for
 reverting
  back earlier. Thanks once again.
 
  With warmest regards
 
  Mike
 
 
  --- On Wed, 25/8/10, Petr PIKAL petr.pi...@precheza.cz wrote:
 
  From: Petr PIKAL petr.pi...@precheza.cz
  Subject: Odp: [R] Finding
   pairs
  To: Mike Rhodes mike_simpso...@yahoo.co.uk
  Cc: r-help@r-project.org
  Date: Wednesday, 25 August, 2010, 6:39
 
  Hi
 
  without other details it is probably impossible to give you any
 reasonable
  advice. Do you have your data already in R? What is their form? Are they

  in 2 columns in data frame? How did you get 

[R] package MuMI

2010-08-25 Thread Marino Taussig De Bodonia, Agnese
Hello,

I am using the package MuMI to run all the possible combinations of variables 
in my full model, and select my best models. When I enter my variables in the 
original model I write them like this

lm(y~ a +b +c +a:b)

However,  MuMI will also use the variable b:a, which I do not want in my 
model.

How do I stop that from happening?

Thank you,

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

2010-08-25 Thread Kay Cecil Cichini

exactly -
thanks a lot, richard!

kay




Zitat von RICHARD M. HEIBERGER r...@temple.edu:


Kay,

doe this do what you want?


dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(x = list(rot = 90, tck=c(1,0))),
ylab=c(y1, y2),
xlab=c(Site 1, Site 2),
strip=FALSE)


On Wed, Aug 25, 2010 at 11:04 AM, Kay Cichini kay.cich...@uibk.ac.atwrote:



hello,

i want to stack two lattice plots beneath each other using one x-axis and
sharing the same text-panels,
like:
#
library(lattice)

y1 - rnorm(100,100,10)
y2 - rnorm(100,10,1)
facs-expand.grid(Sites=rep(c(Site I,Site II),25),Treatment=c(A,B))
pl1-dotplot(y1 ~ facs$Treatment|facs$Sites,
scales = list(x = list(rot = 90, tck=c(1,0
pl2-dotplot(y2 ~ facs$Treatment|facs$Sites,
scales = list(x = list(rot = 90, tck=c(1,0

print(pl1, split=c(1,2,1,2), more=TRUE)
print(pl2, split=c(1,1,1,2))
#

but as said, ideally the plots should be stacked with only the lower plot
giving the x-axis annotation
and only the upper plot with text-panels.

thanks a lot,
kay

-

Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck


--
View this message in context:
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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.htmlhttp://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] Repeat the first day data through all the day. Zoo

2010-08-25 Thread Gabor Grothendieck
On Wed, Aug 25, 2010 at 11:18 AM, skan juanp...@gmail.com wrote:

 # duplicated / na.locf   doesn't work
 it says Error in fix.by(by.x, x) : 'by' must specify valid column(s)

 if I use ifelse instead of ifelse.zoo it works but it gives me a non zoo
 vector.
 Myabe is because my zoo version is older.


They all work:


 library(zoo)
 library(chron)
 z - zoo(1:10, chron(0:9/5))

 # aggregate / na.locf
 z.ag - aggregate(z, as.Date, head, 1)
 na.locf(z.ag, xout = time(z))
(01/01/70 00:00:00) (01/01/70 04:48:00) (01/01/70 09:36:00) (01/01/70 14:24:00)
  1   1   1   1
(01/01/70 19:12:00) (01/02/70 00:00:00) (01/02/70 04:48:00) (01/02/70 09:36:00)
  1   6   6   6
(01/02/70 14:24:00) (01/02/70 19:12:00)
  6   6

 # duplicated / na.locf
 z.na - ifelse.zoo(!duplicated(as.Date(time(z))), z, NA)
 na.locf(z.na)
(01/01/70 00:00:00) (01/01/70 04:48:00) (01/01/70 09:36:00) (01/01/70 14:24:00)
  1   1   1   1
(01/01/70 19:12:00) (01/02/70 00:00:00) (01/02/70 04:48:00) (01/02/70 09:36:00)
  1   6   6   6
(01/02/70 14:24:00) (01/02/70 19:12:00)
  6   6

 # ave - as before
 zz - z
 zz[] - ave(coredata(z), as.Date(time(z)), FUN = function(x) head(x, 1))
 zz
(01/01/70 00:00:00) (01/01/70 04:48:00) (01/01/70 09:36:00) (01/01/70 14:24:00)
  1   1   1   1
(01/01/70 19:12:00) (01/02/70 00:00:00) (01/02/70 04:48:00) (01/02/70 09:36:00)
  1   6   6   6
(01/02/70 14:24:00) (01/02/70 19:12:00)
  6   6

 packageDescription(zoo)$Version
[1] 1.6-4
 R.version.string
[1] R version 2.11.1 Patched (2010-05-31 r52167)

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

2010-08-25 Thread Ben Bolker
 [Sending both to the maintainer and to R-help, in case anyone else has
answers ...]

  I've looked in odfWeave documentation, vignette, and poked around on
the web some, and haven't found answers yet.

  1a. am I right in believing that odfWeave does not respect the
'keep.source' option?  Am I missing something obvious?

  1b. is there a way to set global options analogous to \SweaveOpts{}
directives in Sweave? (I looked at odfWeaveControl, it doesn't seem to
do it.)

  2. I tried to write a Makefile directive to process files from the
command line:

%.odt: %_in.odt
$(RSCRIPT) -e library(odfWeave);
odfWeave(\$*_in.odt\,\$*.odt\);

  This works, *but* the resulting output file gives a warning (The file
'odftest2.odt' is corrupt and therefore cannot be opened.
OpenOffice.org can try to repair the file ...).  Based on looking at
the contents, it seems that a spurious/unnecessary 'Rplots.pdf' file is getting
created and zipped in with the rest of the archive; when I unzip, delete
the Rplots.pdf file and re-zip, the ODT file opens without a warning.
Obviously I could post-process but it would be nicer to find a
workaround within R ...

  3. I find the requirement that all file paths be specified as absolute
rather than relative paths somewhat annoying -- I understand the reason,
but it goes against one practice that I try to encourage for
reproducibility, which is *not* to use absolute file paths -- when
moving a same set of data and analysis files across computers, it's hard
to enforce them all ending up in the same absolute location, which then
means that the recipient has to edit the ODT file.  It would be nice if
there were hooks for read.table() and load() as there are for plotting
and package/namespace loading -- then one could just copy them into the
working directory on the fly.
   has anyone experienced this/thought of any workarounds?
  (I guess one solution is to zip any necessary source files into the archive 
beforehand, 
as illustrated in the vignette.)

  My test files are posted at
http://www.math.mcmaster.ca/~bolker/junk/odftest_in.odt and
http://www.math.mcmaster.ca/~bolker/junk/odftest.odt

 thanks,
Ben Bolker


  sessionInfo()
   
R version 2.11.1 (2010-05-31)
i486-pc-linux-gnu

locale:
 [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=C  LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] odfWeave_0.7.16 XML_3.1-1   lattice_0.19-9

loaded via a namespace (and not attached):
[1] grid_2.11.1

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


Re: [R] frequency, count rows, data for heat map

2010-08-25 Thread Dennis Murphy
Hi:

Here are a couple of ways to render a basic 2D table. Let's call your input
data frame dat:

 names(dat) - c('samp', 'sequen')
 ssTab - as.data.frame(with(dat, table(samp, sequen)))
 ssTab   # data frame version
  samp sequen Freq
1  111abc1
2 1079abc1
3 5576abc1
4  111sdf1
5 1079sdf0
6 5576sdf2
7  111xyz1
8 1079xyz2
9 5576xyz0
 with(dat, table(samp, sequen))   # table version
  sequen
samp   abc sdf xyz
  1111   1   1
  1079   1   0   2
  5576   1   2   0

HTH,
Dennis

On Wed, Aug 25, 2010 at 7:53 AM, rtsweeney tripswee...@gmail.com wrote:


 Hi all,
 I have read posts of heat map creation but I am one step prior --
 Here is what I am trying to do and wonder if you have any tips?
 We are trying to map sequence reads from tumors to viral genomes.

 Example input file :
 111 abc
 111 sdf
 111 xyz
 1079   abc
 1079   xyz
 1079   xyz
 5576   abc
 5576   sdf
 5576   sdf

 How may xyz's are there for 1079 and 111? How many abc's, etc?
 How many times did reads from sample (1079) align to virus xyz.
 In some cases there are thousands per virus in a give sample, sometimes
 one.
 The original file (two columns by tens of thousands of rows; 20 MB) is
 text file (tab delimited).

 Output file:
 abc  sdf  xyz
 111 1  1 1
 1079   1  0 2
 5576   1  2 0

 Or, other ways to generate this data so I can then use it for heat map
 creation?

 Thanks for any help you may have,

 rtsweeney
 palo alto, ca
 --
 View this message in context:
 http://r.789695.n4.nabble.com/frequency-count-rows-data-for-heat-map-tp2338363p2338363.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.


[[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] Change value of a slot of an S4 object within a method.

2010-08-25 Thread Joris Meys
Hi Steve,

thanks for the tip.  I'll definitely take a closer look at your
solution for implementation for future use.  But right now I don't
have the time to start rewriting my class definitions.

Luckily, I found where exactly things were going wrong. After reading
into the documentation about the evaluation in R, I figured out I have
to specify the environment where substitute should look explicitly as
parent.frame(1). I still don't understand completely why exactly, but
it does the job.

Thus :
eval(eval(substitute(expression(obj...@extra[[name]] - value

should become :

eval(
   eval(
  substitute(
 expression(obj...@extra[[name]] - value)
  ,env=parent.frame(1) )
   )
)

Tried it out and it works.

Cheers
Joris

On Wed, Aug 25, 2010 at 6:21 PM, Steve Lianoglou
mailinglist.honey...@gmail.com wrote:
 Hi Joris,

 On Wed, Aug 25, 2010 at 11:56 AM, Joris Meys jorism...@gmail.com wrote:
 Dear all,

 I have an S4 class with a slot extra which is a list. I want to be
 able to add an element called name to that list, containing the
 object value (which can be a vector, a dataframe or another S4
 object)

 Obviously

 setMethod(add.extra,signature=c(PM10Data,character,vector),
  function(object,name,value){
             obj...@extra[[name]] - value
  }
 )

 Contrary to what I would expect, the line :
 eval(eval(substitute(expression(obj...@extra[[name]] - value

 gives the error :
 Error in obj...@extra[[name]] - value : object 'object' not found

 Substitute apparently doesn't work any more in S4 methods...

  I found a work-around by calling the initializer every time, but this
 influences the performance. Returning the object is also not an
 option, as I'd have to remember to assign that one each time to the
 original name.

 Basically I'm trying to do some call by reference with S4, but don't
 see how I should do that. How would I be able to tackle this problem
 in an efficient and elegant way?

 In lots of my own S4 classes I define a slot called .cache which is
 an environment for this exact purpose.

 Using this solution for your scenario might look something like this:

 setMethod(add.extra,signature=c(PM10Data,character,vector),
 function(object, name, value) {
  obj...@.cache$extra[[name]] - value
 })

 I'm not sure what your entire problem looks like, but to get your
 extra list, or a value form it, you could:

 setMethod(extra, signature=PM10Data,
 function(object, name=NULL) {
  if (!is.null(name)) {
    obj...@.cache$extra[[name]]
  } else {
    obj...@.cache$extra
 })

 ... or something like that.

 The last thing you have to be careful of is that you nave to make sure
 that each new(PM10Data) object you have initializes its *own* cache:

 setClass(PM10Data, representation(..., .cache='environment'))
 setMethod(initialize, PM10Data,
 function(.Object, ..., .cache=new.env()) {
  callNextMethod(.Object, .cache=.cache, ...)
 })

 Does that make sense?

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


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


Re: [R] How to obtain seed after generating random number?

2010-08-25 Thread Greg Snow
If you find yourself doing things like this often, but don't want to explicitly 
set the seed, or save seeds before simulating, then you can run the following 
code (or put it into .Rprofile or similar):

.Last.Random.seed - .Random.seed

addTaskCallback( function(expr, val, ok, visible){
if(!isTRUE( all.equal(.Last.Random.seed, .Random.seed)) ) {
.Last.Random.seed - .Random.seed
}
TRUE
})

Then the previous seed will be stored in .Last.Random.seed and you can restore 
the seed to be able to rerun the same values, e.g.:

 rnorm(10)
 [1] -0.28361138  0.86951931 -0.54435528  0.62880324 -1.42233446 -1.22751263
 [7] -1.67410552  0.08439848 -0.20612566  1.44187164
 .Random.seed - .Last.Random.seed
 rnorm(10)
 [1] -0.28361138  0.86951931 -0.54435528  0.62880324 -1.42233446 -1.22751263
 [7] -1.67410552  0.08439848 -0.20612566  1.44187164
 rnorm(10)
 [1] -0.0417821  1.3537545  1.9452253 -0.4909382  0.3884391 -0.8448933
 [7]  0.7379904 -1.0797603 -1.0264739  0.2887934

The above code only keeps the most recent seed, but the code could be modified 
to store a longer history if that were desired.


If you want to set your own seeds (a little easier to save, pass to others), 
but find integers to unimaginative, then look at the char2seed function in the 
TeachingDemos package.

Hope this helps,

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


 -Original Message-
 From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-
 project.org] On Behalf Of Bogaso Christofer
 Sent: Tuesday, August 24, 2010 11:12 AM
 To: r-help@r-project.org
 Subject: [R] How to obtain seed after generating random number?
 
 Dear all, I was doing an experiment to disprove some theory therefore
 performing lot of random simulation. Goal is to show the audience that
 although something has very rare chance to occur but it doesn't mean
 that
 event would be impossible.
 
 
 
 In this case after getting that rare event I need to show that same
 scenario
 for multiple times to explain other audience. Hence I need to somehow
 save
 that seed which generates that random numbers after doing the
 experiment.
 However as it is very rare event it is not very practical to start with
 a
 fixed seed and then generate random numbers. Hence I am looking for
 some way
 which will tell me about that corresponding seed which was responsible
 to
 generate that particular series of random numbers responsible for
 occurrence
 of that rare event.
 
 
 
 In short, I need to know the seed ***after*** generating the random
 numbers.
 
 
 
 Is there any possibility to know this?
 
 
 
 Thanks and regards,
 
 
   [[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] nls self starting function

2010-08-25 Thread Marlin Keith Cox
I need the simple function for the following set of data.  This is a toy
version of my data, but the error is persistent in both.  To compare with
excel, I would  just 1) format trendline, 2) display equation and
R-squared on chart.

I obviously tried to use a nls and wanted to use the self startup function
SSasymp

rm(list=ls())
t - c( 0,.1,.2,.7,.4,.5,.6,1,.8,.4,1.2,1.3,1.4,1.5,1.6,
1.6,1.7,1.8,1.9,2.0, 3,6,11)
y-c(2,3.5,3.01,3.09,3,3.25,3.36,3.49,3.64, 3.81, 4.44, 4.69, 4.96, 5.25, 6,
6.1, 5.89, 6.24, 6.61, 7.00,
12.00, 39, 124.00)
plot(y~t ,xlab=t, ylab=y)
model-nls(y~SSasymp(t,a,b))
summary(model)
Thanks in advance,
keith


-- 
M. Keith Cox, Ph.D.
Alaska NOAA Fisheries, National Marine Fisheries Service
Auke Bay Laboratories
17109 Pt. Lena Loop Rd.
Juneau, AK 99801
keith@noaa.gov
marlink...@gmail.com
U.S. (907) 789-6603

[[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] lattice help required

2010-08-25 Thread Kay Cichini

... i added relation=free to account for diffferent ranges of y1 and y2:

dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
 outer=TRUE,
 scales = list(y = list(relation=free), x = list(rot = 90,
tck=c(1,0))),
 ylab=c(y1, y2),
 xlab=c(Site 1, Site 2),
 strip=FALSE)  

but then i get four different y-lims.
more suited there should be only two - one for each response.
can someone tell me how i have to change the code?

best,
kay



 Kay,

 doe this do what you want?


 dotplot(y1+y2 ~ facs$Treatment|facs$Sites,
 outer=TRUE,
 scales = list(x = list(rot = 90, tck=c(1,0))),
 ylab=c(y1, y2),
 xlab=c(Site 1, Site 2),
 strip=FALSE)


 On Wed, Aug 25, 2010 at 11:04 AM, Kay Cichini
 kay.cich...@uibk.ac.atwrote:


 hello,

 i want to stack two lattice plots beneath each other using one x-axis and
 sharing the same text-panels,
 like:
 #
 library(lattice)

 y1 - rnorm(100,100,10)
 y2 - rnorm(100,10,1)
 facs-expand.grid(Sites=rep(c(Site I,Site
 II),25),Treatment=c(A,B))
 pl1-dotplot(y1 ~ facs$Treatment|facs$Sites,
 scales = list(x = list(rot = 90, tck=c(1,0
 pl2-dotplot(y2 ~ facs$Treatment|facs$Sites,
 scales = list(x = list(rot = 90, tck=c(1,0

 print(pl1, split=c(1,2,1,2), more=TRUE)
 print(pl2, split=c(1,1,1,2))
 #

 but as said, ideally the plots should be stacked with only the lower plot
 giving the x-axis annotation
 and only the upper plot with text-panels.

 thanks a lot,
 kay

 -
 
 Kay Cichini
 Postgraduate student
 Institute of Botany
 Univ. of Innsbruck
 

 --
 View this message in context:
 http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338382.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.htmlhttp://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.




-

Kay Cichini
Postgraduate student
Institute of Botany
Univ. of Innsbruck


-- 
View this message in context: 
http://r.789695.n4.nabble.com/lattice-help-required-tp2338382p2338641.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] lattice help required

2010-08-25 Thread RICHARD M. HEIBERGER
The multiple y axes are protecting you in this situation.


z - cbind(rnorm(100,c(1,10),1), rnorm(100,c(20,30),1))
dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(
  y = list(
relation=free)),
ylab=c(y1, y2),
xlab=c(Site 1, Site 2),
strip=FALSE,
main=problem)

dotplot(z[,1]+z[,2] ~ facs$Treatment|facs$Sites,
outer=TRUE,
scales = list(
  y = list(
relation=free,
limits=list(c(-5,13),c(-5,13),c(18,32),c(18,32,
ylab=c(y1, y2),
xlab=c(Site 1, Site 2),
strip=FALSE, main=protected)

For more control (such as suppressing the y-tick labels in the right-hand
column,
I recommend Deepayan Sarkar's book.

Rich

[[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] Secant Method Convergence (Method to replicate Excel XIRR/IRR)

2010-08-25 Thread Adrian Ng
Hi,

I am new to R, and as a first exercise, I decided to try to implement an XIRR 
function using the secant method.  I did a quick search and saw another posting 
that used the Bisection method but wanted to see if it was possible using the 
secant method.

I would input a Cash Flow and Date vector as well as an initial guess.  I 
hardcoded today's initial date so I could do checks in Excel.  This code seems 
to only converge when my initial guess is very close to the correct IRR.

Maybe I have some basic errors in my coding/logic? Any help would be greatly 
appreciated.

The Wikipedia article to secant method and IRR: 
http://en.wikipedia.org/wiki/Internal_rate_of_return#Numerical_solution

Thanks!



ANXIRR - function (cashFlow, cfDate, guess){
cfDate-as.Date(cfDate,format=%m/%d/%Y)
irrprev - c(0); irr- guess


pvPrev- sum(cashFlow)
pv- 
sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,2010-08-24,units=days))/360)))
print(pv)
print(Hi)


while (abs(pv) = 0.001) {
t-irrprev; irrprev- irr;
irr-irr-((irr-t)*pv/(pv-pvPrev));
pvPrev-pv;

pv-sum(cashFlow/((1+irr)^(as.numeric(difftime(cfDate,2010-08-24,units=days))/365)))
print(irr);print(pv)
}
}





Please consider the environment before printing this e-mail.

[[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] Merging two data set in R,

2010-08-25 Thread Danne Mora
try the following merge command.

merge(A,B, by = intersect(names(A), names(B)), all.x=FALSE, all.y=FALSE)
or
merge(A,B, by = ID, all.x=FALSE, all.y=FALSE)

Dannemora

On Wed, Aug 25, 2010 at 5:35 AM, Mangalani Peter Makananisa 
pmakanan...@sars.gov.za wrote:

 Dear R Gurus,



 I am currently working on the two dataset ( A and B), they both have the
 same fields:ID , REGION, OFFICE, CSTART, CEND, NCYCLE, STATUS and
 CB.

 I want to merge the two data set by ID. The problem I have is that the
 in data A, the ID's are unique. However in the data set B, the ID's are
 not unique, thus some repeat themselves.



 How do I the merge or retrieve the common ones?

 Please advise.



 Kind Regards



 Peter



 South Africa

 +27 12 422 7357

 +27 82 456 4669
















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

[[alternative HTML version deleted]]

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


[[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] Showing a line function on the ploat area

2010-08-25 Thread Neta

I have a collection of results. I use R to get the linearization presented.
how can I get R to show the equation and R^2 value on the plot area next to
the graph?
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Showing-a-line-function-on-the-ploat-area-tp2338389p2338389.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] Plot bar lines like excel

2010-08-25 Thread abotaha

Great ..
thanks for the to much help and i too appreciate hel p and explanation 

Cheers




-- 
View this message in context: 
http://r.789695.n4.nabble.com/Plot-bar-lines-like-excel-tp2337089p2338158.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] GLM outputs in condensed versus expanded table

2010-08-25 Thread francogrex

Hi I'm having different outputs from GLM when using a condensed table
V1  V2  V3  Present Absent
0   0   0   3   12
0   0   1   0   0
0   1   0   0   0
0   1   1   1   0
1   0   0   7   20
1   0   1   0   0
1   1   0   3   0
1   1   1   6   0


resp=cbind(Present, Absent)
glm(resp~V1+V2+V3+I(V1*V2*V3),family=binomial)
Deviance Residuals: 
[1]  0  0  0  0  0  0  0  0
 etc and also coefficients...

And when using the same but expanded table

V1  V2  V3  condition (1 present 0 abscent)
Id1 1   0   0   1
id2 1   1   1   1
... etc
glm(condition~V1+V2+V3+I(V1*V2*V3),family=binomial)
Deviance Residuals: 
Min  1Q  Median  3Q Max  
  -0.7747317  -0.7747317  -0.6680472   0.0001315   1.7941226 
and also coefficients are different from above.

What could I be doing wrong?


-- 
View this message in context: 
http://r.789695.n4.nabble.com/GLM-outputs-in-condensed-versus-expanded-table-tp2338407p2338407.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] NewJerseyR Meeting

2010-08-25 Thread Sarah Lewis
Mango Solutions are proud to announce the first NewJerseyR Meeting  to be held 
on the 16th September 2010.

Please see details below:

NewJerseyR

Date:   Tuesday 16th September
Time:   6.30pm - 9.30pm
Venue:  Renaissance Woodbridge Hotel - Iselin, New Jersey

The agenda has yet to be finalised, I shall send details out as soon as they 
have been confirmed.

Please contact: newjers...@mango-solutions.com  to be added to the NewJerseyR 
mailing list to receive meeting updates and R news.

To register for this event please contact: newjers...@mango-solutions.com  

We currently hold 2 very successful R meetings - LondonR and BaselR. NewJerseyR 
will be our first R meeting in America. NewJerseyR will be a free informal 
event, held quarterly and intended to serve as a platform for all local (and 
regional) R users to present and exchange their experiences and ideas around 
the usage of R. A typical meeting will consist of 3-4 talks of about 20-25 min 
to give plenty of room for sharing your R experiences, discussions and exchange 
of ideas. 

For presentations from our other R meetings can be found at:
www.londonr.org   and  www.baselr.org  
The NewJerseyR website will be up and running shortly.

For further information, please contact newjers...@mango-solutions.com 


Sarah Lewis

Hadley Wickham, Creator of ggplot2 - first time teaching in the UK. 1st - 2nd  
November 2010. 
To book your seat please go to http://mango-solutions.com/news.html 
T: +44 (0)1249 767700 Ext: 200
F: +44 (0)1249 767707
M: +44 (0)7746 224226
www.mango-solutions.com 
Unit 2 Greenways Business Park 
Bellinger Close
Chippenham
Wilts
SN15 1BN
UK 


LEGAL NOTICE
This message is intended for the use o...{{dropped:9}}

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


Re: [R] how to plot y-axis on the right of x-axis

2010-08-25 Thread Antonio Olinto

Dear Elaine

I'm developing a code to make a 3 y-axes plot. It may can help you.

Also your question leads to mine, posted yesterday: How does R  
calculate the interval between tick-marks.


Below follows the code I'm developing.

Data:

ANO CKG NUP NDE
200526352158150025014
200632514789210035024
200727458126180030254
200828568971150043254
200932564789300060320

Code:

# I use comma as decimal symbol so I use dec=,
dat.caupde - read.delim(clipboard,dec=,,header=T)

attach(dat.caupde)
summary(dat.caupde)
  ANOCKGNUPNDE
 Min.   :2005   Min.   :26352158   Min.   :1500   Min.   :25014
 1st Qu.:2006   1st Qu.:27458126   1st Qu.:1500   1st Qu.:30254
 Median :2007   Median :28568971   Median :1800   Median :35024
 Mean   :2007   Mean   :29491767   Mean   :1980   Mean   :38773
 3rd Qu.:2008   3rd Qu.:32514789   3rd Qu.:2100   3rd Qu.:43254
 Max.   :2009   Max.   :32564789   Max.   :3000   Max.   :60320

# here I indicate the limits of each axes and calculate their ranges
y1min - 26
y1max - 33
y1range - y1max-y1min
y2min - 1500
y2max - 3000
y2range - y2max-y2min
y3min - 25
y3max - 61
y3range - y3max-y3min

# making the plot
# y1range*((NUP-y2min)/y2range)+y1min calculates the proportion to  
between axes


par(mar=c(6, 6, 2,12))
plot(CKG/100, ylim=c(y1min,y1max), ylab=landings (1000 t),  
type=l,col=red,las=0, cex.axis=1.2,cex.lab=1.4,xaxt=n,xlab=)

points( y1range*((NUP-y2min)/y2range)+y1min,type=l,col=blue)
points( y1range*((NDE/1000-y3min)/y3range)+y1min,type=l,col=darkgreen)

# the number 1 in axis(4,at=seq(y1min,y1max,1 ... is the interval  
between each tick-mark

axis(1,at=1:5,labels=as.character(ANO),las=2,cex.axis=1.2)
axis(4,at=seq(y1min,y1max,1),labels=as.integer(seq(y2min,y2max,y2range/y1range)),las=0,cex.axis=1.2)
axis(4,at=seq(y1min,y1max,1),  
labels=as.integer(seq(y3min,y3max,y3range/y1range)),cex.axis=1.2,  
las=0,line=5)

mtext(nº of fishing boats,4,3,cex=1.4)
mtext(nº of fishing trips (x1000),4,8, cex=1.4)
legend(4,33,c(L,Nº FB,Nº  
FT),bty=n,lty=c(1,1,1),col=c(red,blue,darkgreen),cex=1.4)


Now I want to replace the number 1 by the formula used to calculate  
the interval between tick-marks. Why, with the range 26-33, R choose  
unitary intervals for the y axis (26, 27, 28 ...)?


All the best,

Antonio Olinto


Citando elaine kuo elaine.kuo...@gmail.com:


Dear List,

I have a richness data distributing across 20 N to 20 S latitude. (120 E-140
E longitude).


I would like to draw the richness in the north hemisphere and a regression
line in the plot
(x-axis: latitude, y-axis: richness in the north hemisphere).
The above demand is done using plot.

Then, south hemisphere richness and regression are required to be generated
using
the same y-axis above but an x-axis on the left side of the y-axis.
(The higher latitude in the south hemisphere, the left it would move)

Please kindly share how to design the south plot and regression line for
richness.
Also, please advise if any more info is in need.

Elaine

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






Webmail - iBCMG Internet
http://www.ibcmg.com.br

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/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 remove all objects except a few specified objects?

2010-08-25 Thread Cheng Peng

Thanks to De-Jian and Peter. Peter's way is neat and cool! 
-- 
View this message in context: 
http://r.789695.n4.nabble.com/How-to-remove-all-objects-except-a-few-specified-objects-tp2335651p2338221.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] Problem with clusterCall, Error in checkForRemoteErrors(lapply(cl, recvResult)) :

2010-08-25 Thread telm8

Hi all,

I am trying to use snow package to do a parallel MCMC. I have read a few
guides and articles, the following is that I came up with.

When I run it I got the error message:
Error in checkForRemoteErrors(lapply(cl, recvResult)) : 
  4 nodes produced errors; first error: could not find function ui.Next

The data is a longitudinal data with few repeated readings on a number of
individuals. However, the response is organised in a vector rather than a
matrix. E.g. (y_11 , y_12 , y_13 , y_14 , y_21 , y_22 , ... , y_n1 , ... ,
y_nTn )^T

X , Z are covariates.

beta is a matrix of coefficients associated with X

C is the latent class membership indicator

sigma.sq is the diagonal elements of the covariance matrix of u, which is a
random effect parameter and the one I am trying to sample with the function.

sd is the diagonal elements of the covariance matrix of the proposal
distribution.

The following is the particular function:

ui.Full.Sample.Cluster = function( levels , Y , X , beta , Z , C , sigma.sq
, sd , burnin , iteration )
{
cluster = makeCluster( 4 , type = MPI);
clusterSetupRNG( cluster );

patients = unique( levels );

q = length( X[1,] );

u = ui.Ini( q , length( patients ) );

n = levels[ length(patients) ];

marker = levels == patients[1];

y = Y[ marker ];
x = X[ marker, ];
z = Z[ marker, ];

u[1,] = ui.1.Sample( u[1,] , y , x , beta , z , C[1] , sigma.sq , sd ,
burnin , iteration )$uFinal;

for( i in 2:n)
{
marker = levels == patients[i];

y = Y[ marker ];
x = X[ marker, ];
z = Z[ marker, ];

print( i );

u[i, ] = clusterCall( cluster , ui.1.Sample, u[i,] , y , x , 
beta , z ,
C[i] , sigma.sq , sd , burnin , iteration )$uFinal;

print( i );
}

stopCluster( cluster );

u;
}



If anyone could help that would be much appreciated! But big thanks for
taking your time to read the post!

TelM8
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Problem-with-clusterCall-Error-in-checkForRemoteErrors-lapply-cl-recvResult-tp2338375p2338375.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] What does this warning message (from optim function) mean?

2010-08-25 Thread Sally Luo
Hi R users,
I am trying to use the optim function to maximize a likelihood funciton, and
I got the following warning messages.
Could anyone explain to me what messege 31 means exactly?  Is it a cause for
concern?
Since the value of convergence turns out to be zero, it means that the
converging is successful, right?
So can I assume that the parameter estimates generated thereafter are
reliable MLE estimates?
Thanks a lot for your help.

Maomao

 p-optim(c(0,0,0), f, method =BFGS, hessian =T, y=y,X=X,W=W)

There were 31 warnings (use warnings() to see them)

 warnings()

Warning messages:

1: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

2: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

3: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

4: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

5: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

6: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

7: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

8: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

9: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

10: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

11: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

12: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

13: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

14: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

15: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

16: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

17: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

18: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

19: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

20: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

21: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

22: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

23: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

24: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

25: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

26: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

27: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

28: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

29: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

30: In log(det(I_N - pd * wd - po * wo - pw * ww)) : NaNs produced

31: In if (hessian) { ... :

  the condition has length  1 and only the first element will be used

 p$counts

function gradient

 148   17

 p$convergence

[1] 0

[[alternative HTML version deleted]]

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


[R] pairwise correlation of many samples

2010-08-25 Thread Moon Eunyoung
Dear list,



I have a csv file like below, and I have similar 150 files (responses from
150 samples) in a folder.

I’d like to get the mean score of pairwise correlation score among 150
respondents by hour(0~20hour).

All I can think of is bring up two files and get the pairwise correlation
score and repeating this (e.g. (p1,p2),(p1,P3)…(P4,P5)..)and get the total
mean of correlation score.

However, since there’re 150 samples I can’t think of doing this manually..







Hour

Response



0

-0.11703



0.016667

-0.10427



0.03

-0.091347



0.05

-0.078313



0.07

-0.065226



0.08

-0.052147



0.1

-0.039136



….

….





Can anyone please point out the best way to do? .. It’s been only a month
since I started to use R



Thank you!

[[alternative HTML version deleted]]

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


Re: [R] Problem with clusterCall, Error in checkForRemoteErrors(lapply(cl, recvResult)) :

2010-08-25 Thread telm8

Oh, I forgot to say, ui.Next() is a function I defined to sample from the
proposal distribution given the current state.
-- 
View this message in context: 
http://r.789695.n4.nabble.com/Problem-with-clusterCall-Error-in-checkForRemoteErrors-lapply-cl-recvResult-tp2338375p2338378.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] multivariate analysis of variance

2010-08-25 Thread celal arslan
Hi,

can anyone tell me how i may find out the between-group variance and 
within-group variance for multivariate case?

I have 6 Groups and 73 Variables. (with MANOVA ? wie) 
dim(data)

[1] 2034   76


Thanks

Celal 




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