Re: [R] return one vector that is the sum of a list of vectors

2012-04-09 Thread David Winsemius


On Apr 9, 2012, at 12:01 AM, Frank Tamborello wrote:


Dear R Help,

I am attempting to write a function that  takes a list of variable  
groups and a vector of numbers (e.g., jtsv would indicate one  
group of variables, jtsv1, jtsv2, etc, each of which name a  
variable in a data frame), and returns a list of variable group  
sums. For example, given list(jtsv, ptsv) and c(1:10), where  
jtsv1, etc, are all numeric vectors of the same length, the function  
should return a list of jtsv and ptsv, where
jtsv - jtsv1 + jtsv2 + jtsv3 + jtsv4 + jtsv5 + jtsv6 + jtsv7 +  
jtsv8 + jtsv9 + jtsv10
ptsv - ptsv1 + ptsv2 + ptsv3 + ptsv4 + ptsv5 + ptsv6 + ptsv7 +  
ptsv8 + ptsv9 + ptsv10


So far I've used a for loop to paste together the names of the  
variables as character vectors, so that I can at least name the  
variables that I want to get, but I'm stuck on how to add together  
the numeric vectors that those character vectors are ultimately  
supposed to name. It seems like I should be able to map the  
primitive + operator onto a list or vector that contains the  
variable names. If that's a good way to go, what would that look  
like? Or if that's not a good way to go, would someone please point  
me in a good direction?




Perhaps ( depending on detail not provided):

jtsv - rowSums( dfrm[ , grep(jtsv, names(dfrm))] )

--

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] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread Guaramy
Sorry but i didn't understand.

Thanks for your answer.

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-integrate-int-fn-lower-0-upper-Inf-evaluation-of-function-gave-a-result-of-wrong-length-tp4541250p4542101.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] Creating a dummy variable

2012-04-09 Thread arunkumar1111
HI

I want to create a dummy variable for a dataset which has various dates Eg
I've one year data. The user can choose any date range

startdate : 2012-01-01
enddate: 2012-02-01

startdate: 2012-06-01
enddate=2012-07-01

the columns of dataset= X ,Y ,date
   

Thanks and re

-
Thanks in Advance
Arun
--
View this message in context: 
http://r.789695.n4.nabble.com/Creating-a-dummy-variable-tp4542150p4542150.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] Summary Statistics Help

2012-04-09 Thread bobo
I've got this solved via Talks Stat

mod.1-lm(Patents~FHouse, data=datpat) 
summary(mod.1) 
anova(mod.1)
xtable(mod.1)

--
View this message in context: 
http://r.789695.n4.nabble.com/Summary-Statistics-Help-tp4541923p4542103.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] Problems with CRAN package Ryacas

2012-04-09 Thread Jeff Newmiller
check your firewall?
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

Kjetil Halvorsen kjetilbrinchmannhalvor...@gmail.com wrote:

First:
 sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-pc-linux-gnu (64-bit)

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=en_US.UTF-8LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=C 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] Ryacas_0.2-11 XML_3.9-4 pnmath_0.0-3  MASS_7.3-17

loaded via a namespace (and not attached):
[1] compiler_2.15.0 fortunes_1.5-0  tools_2.15.0

Operating system is debian (wheezy+some sid)

I have installed yacas version 1.3.1 via synaptic.

Then:

library(Ryacas)
 yacas(expression(Factor(x^2-1))) # copied from yacas help page
[1] Starting Yacas!
Accepting requests from port 9734
Error in socketConnection(host = 127.0.0.1, port = 9734, server =
FALSE,  :
  cannot open the connection
In addition: Warning message:
In socketConnection(host = 127.0.0.1, port = 9734, server = FALSE,  :
  127.0.0.1:9734 cannot be opened



¿Any ideas?

Kjetil

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

2012-04-09 Thread Frank Schwidom
On Thu, Sep 22, 2011 at 06:46:23PM +, Joseph Boyer wrote:
 Would someone be so kind as to provide example code where they use the 
 suggestions argument in  the rgba function
 In genalg? I can't get it to work.
 
 The following code works just fine:
 
 GenFit -rbga(Lower, Upper, evalFunc = evaluate)
 
 Lower and Upper are each numeric vectors with 7 elements. Evaluate is an 
 objective function.
 However, when I want to use a suggested chromosome, I get an error message. 
 My code is
 
 start - c(1,0.1,10, 100,1,100,1)
 
 suggestions - list(start)
 
 GenFit -rbga(Lower, Upper, suggestions = suggestions, evalFunc = evaluate)
 
 The error message is:
 
 Error in 1:suggestionCount : argument of length 0
 
 Thanks.

rbga( c( -1, -1), c( 1, 1), evalFunc= function( string){ sum(
string**2)**0.5}, suggestions=t( as.data.frame( c( 0, 0

Regards

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

2012-04-09 Thread ali_protocol
I am interested in the difference of 2 data:

mat1= c(2.2, 2.3, 2.2,2.5)
mat2= c(2.6, 2.8, 2.7,2.4)

mat= mat2-mat1

I perform an action on both mat1 and mat2, and I get
mat1prime and mat2prime:

mat1prime= c(2.5, 2.5, 2.3,2.5)
mat2prime= c(2.6, 2.8, 2.7,2.6)

matprime= mat2prime-mat1prime


I want to check wether mat and matprime have the same mean.
How should I do this in R?

--
View this message in context: 
http://r.789695.n4.nabble.com/Comparing-2-means-pls-help-tp4542237p4542237.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] slanted stacked bar graphs?

2012-04-09 Thread Barry Rowlingson
On Mon, Apr 9, 2012 at 7:29 AM, Susanna Makela
susanna.m.mak...@gmail.com wrote:
 Hello R users,

 I would like to generate slanted stacked bar graphs like those on
 the bottom of pages 1 and 2 in this document:
 http://www.wssinfo.org/fileadmin/user_upload/resources/JMP-Snapshot-SWA-HLM.pdf
 . I've also attached the file to this email (pdf). Does anyone know if
 this is possible in R? I have tried googling and searching the R help
 archives, and it seems like ggplot2 might be able to make such graphs,
 but I'm not familiar enough with graphics in R to know for sure.

 (I personally don't feel that these slanted bar graphs - not sure if
 they have an actual name - convey the intended information very well,
 but I have to try and make them all the same. However, I am open to
 alternative suggestions for visualizing similar data if anyone has
 ideas.)


These exact charts have been critiqued on the Junk Charts blog:

http://junkcharts.typepad.com/junk_charts/2010/02/cousin-misfit.html

 and you'll even find some ggplot code in the comments for doing them.
If you still want to...

 I just did a google image search for 'ggplot stacked' and there they were.

Barry

-- 
blog: http://geospaced.blogspot.com/
web: http://www.maths.lancs.ac.uk/~rowlings
web: http://www.rowlingson.com/
twitter: http://twitter.com/geospacedman
pics: http://www.flickr.com/photos/spacedman

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


Re: [R] Drawing a line in xyplot

2012-04-09 Thread John Maindonald
PS: The following shows possibilities that are available using latticeExtra 
layering:

## Best make type and attend factors
xdat = data.frame(mortality =c(5, 
8,7,5,8,10,11,6,4,5,20,25,27,30,35,32,28,21,20,34,11,15,18,12,15,12,10,15,19,20),
 type=
factor(c(1, 1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3)),
attend = factor(c(1, 
0,1,1,1,1,0,0,0,1,0,1,0,0,0,0,1,0,0,0,1,1,1,0,1,0,1,0,0,0)))


gph - xyplot (mortality ~ attend|type, pch=16,
 data = xdat, aspect=2:1,layout=c(3,1))

one4all - layer(avy - median(xdat$mortality),
  panel.segments(0.1, avy, 0.3, avy, col=red,lwd=4),
  panel.segments(0.7, avy, 1, avy, col=red,lwd=4))

medbytype - layer(avy - median(y),
  panel.segments(0.1, avy, 0.3, avy, col=red,lwd=4),
  panel.segments(0.7, avy, 1, avy, col=red,lwd=4))

interact - layer(panel.average(x, y, fun=median, col='red', lwd=4))

Compare

gph + one4all
(shows overall median lines in all 3 panels)

gph + medbytype
(shows separate median lines for the separate panels)

gph+interact
(gives a form of interaction plot)

NB x (if its values are used) and y are local to the individual panel

NB also that layer() accepts a data argument.  The following is an
alternative way to calculate one4all:

one4all - layer(data=xdat, avy - median(mortality),
  panel.segments(0.1, avy, 0.3, avy, col=red,lwd=4),
  panel.segments(0.7, avy, 1, avy, col=red, lwd=4))

John Maindonald.


 This can be simplified by using the layering abilities that Felix Andrews 
 made available in latticeExtra.  These are too little known.  These pretty
 much make it unnecessary to resort to trellis.focus(), at least in such
 cases as this.  These layering abilities are too little known:
 
 library(latticeExtra)
 x11(height=8,width=11)
 xdat = data.frame(mortality =c(5, 
 8,7,5,8,10,11,6,4,5,20,25,27,30,35,32,28,21,20,34,11,15,18,12,15,12,10,15,19,20),
  type=
 c(1, 1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3), attend = 
 c(1, 0,1,1,1,1,0,0,0,1,0,1,0,0,0,0,1,0,0,0,1,1,1,0,1,0,1,0,0,0))
 gph - xyplot ( mortality ~ attend|type,
 panel=function(x,y)
 {
 panel.grid(lty=5)
 panel.xyplot(x,y,pch=16,jitter.x=TRUE,col=1)
 for(i in 1:3)
 {
 panel.segments(x0=c(0.7, 1.7),
 x1=c(1.3, 2.3),
 y0=with(xdat, c(tapply(mortality[type==i], attend[type==i], median))),
 y1=with(xdat, c(tapply(mortality[type==i], attend[type==i],
 median))),col=red,lwd=4)
 } 
 },
 data = xdat, aspect=2:1,layout=c(3,1))
 
 ## Now create a layer object that will add the further segments.
 addlayer - layer(panel.segments(x0=1,y0=20,x1=1.2, y1=20),
 panel.segments(x0=2,y0=30,x1=2.2, y1=30))
 
 gph+addlayer
 
 The code that produces the object gph would also be simpler
 and easier to follow if some relevant part was separated out into 
 a separate layer.
 
 See also my notices on layering of lattice objects at:
 http://www.maths.anu.edu.au/%7Ejohnm/r-book/add-graphics.html
 
 John Maindonald email: john.maindon...@anu.edu.au
 phone : +61 2 (6125)3473fax  : +61 2(6125)5549
 Centre for Mathematics  Its Applications, Room 1194,
 John Dedman Mathematical Sciences Building (Building 27)
 Australian National University, Canberra ACT 0200.
 http://www.maths.anu.edu.au/~johnm
 
 On 08/04/2012, at 8:00 PM, r-help-requ...@r-project.org wrote:
 
 From: David Winsemius dwinsem...@comcast.net
 Subject: Re: [R] Drawing a line in xyplot
 Date: 8 April 2012 2:17:22 PM AEST
 To: wcheckle wchec...@jhsph.edu
 Cc: r-help@r-project.org
 
 
 
 On Apr 7, 2012, at 10:29 PM, wcheckle wrote:
 
 Thank you David, the bwplot option does what I need:
 
 x11(height=8,width=11)
 bwplot ( mortality~ attend|type,
 pch=95,cex=5,col=2,
 par.settings=list(
 box.rectangle = list(col = transparent),
 box.umbrella = list(col = transparent),
 plot.symbol = list(col = transparent)
 ),
 panel=function(x,y,...){
 panel.grid(lty=5)
 panel.xyplot(x,y,pch=16,jitter.x=TRUE,col=1)
 panel.bwplot(x,y,...)
 },
 data = x,aspect=2:1,layout=c(3,1))
 
 
 However, I am interested in also learning how to do it in xyplot as well.  I
 wasn’t able to follow the last two set of instructions (condition on
 packet.number and loop over segments), wondering if I can ask for your help
 for the correct code (my attempt resulted in all three mean lines within
 each panel):
 
 x11(height=8,width=11)
 xyplot ( mortality ~ attend|type,
 panel=function(x,y)
 {
 panel.grid(lty=5)
 panel.xyplot(x,y,pch=16,jitter.x=TRUE,col=1)
 for(i in 1:3)
 {
 panel.segments(x0=c(0.7, 1.7),
 x1=c(1.3, 2.3),
 y0=c(tapply(x$mortality[x$type==i], x$attend[x$type==i], median)),
 y1=c(tapply(x$mortality[x$type==i], x$attend[x$type==i],
 median)),col=red,lwd=4)
 }   
 },
 data = x,aspect=2:1,layout=c(3,1))
 
 thank you again. I also found your info on data.frame useful.
 
 I haven't figured out how to do it with the interaction of 'attend' and 
 ''type' from inside xyplot. I'm thinking I might be able to succeed using 
 trellis.focus() to address separate columns within a 

[R] Stepwise procedure with force.in command

2012-04-09 Thread Hien Nguyen

Dear R-helpers,

I am trying to do a stepwise procedure in which I want to force some 
variables in the model. I have searched around and it seems that only 
leaps package allows to force the variable in the stepwise procedure. I 
use the leaps package and use the regsubsets(lm1, force.in = 1, data) to 
force 1 variable in the model. However, the force.in command only allow 
me to force 1 variable in. I want to force several variables in the 
procedure. How could I force them in?


Thanks a lot,

Hien

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 get the confidence interval of area under the time dependent roc curve

2012-04-09 Thread Terry Therneau

On 04/07/2012 05:00 AM, r-help-requ...@r-project.org wrote:

  It is possible to calculate the c-index for time dependent outcomes (such as 
disease) using the survivalROC package in R. My question is : is it possible to 
produce a p-value for the c-index that is calculated (at a specific point in 
time)?How to get the confidence interval of area under then time dependent roc 
curve (or by bootstrap)?
The current version of survival computes the c-index and it's standard 
error for time-dependent covariates.

Use summary(cox model fit) or the survConcordance function.

Terry T.

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


Re: [R] Order sapply

2012-04-09 Thread MSousa
Good Afternoon,

 I get it, my question was really how the data were organized, I thought
I was doing something wrong.
  
Thanks

--
View this message in context: 
http://r.789695.n4.nabble.com/Order-sapply-tp4537496p4542518.html
Sent from the R help mailing list archive at Nabble.com.

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


[R] Question on harmonic (Fourier) analysis of sinusoidal time series

2012-04-09 Thread holden999
Hello,

I will try to explain the problem, sorry if  it will be a little long...

I'm using R to analyze results of cyclic mechanical testing, like this:
- apply quasi-sinusoidal load
- measure quasi-sinusoidal vertical and horizontal deformations

(quasi-sinusoidal load means that load should be sinusoidal, but testing
machine puts in some noise...)
I enclose a sample of data at the end of this message.

I used harmonic regression to fit the data: 

force = constant + [1st order Fourier harmonic] + [2nd order Fourier
harmonic] + higher order harmonics...  + error =

= v0 + [a1 sin (2*pi*f*t) +  b1 cos (2*pi*f*t) ] + [a2 sin (2*2*pi*f*t) + 
b2 cos (2*2*pi*f*t)] + ... + error

where f is the frequency of the sinusoidal load applied during the test.

What I'm expecting to find is:
- v0 (average value)
- a1, b1 high values of coefficients of the first fundamental harmonic
that the machine is applying 
- an, bn  low values of coefficients of the higher  harmonics 
- some noise ~ NID

I used lm() to fit the linear model using the code suggested in
Introductory Time Series with R
by Paul S.P. Cowpertwait (2006 SPRINGER SCIENCE+BUSINESS MEDIA, LLC.).
Started with 1st harmonic component and then increased up to 6.

Here's the results (6 harmonic components, but results are similar with less
harmonics):

Call:
lm(formula = force.fit ~ s + c)

Residuals:
   Min 1Q Median 3QMax 
-0.0147651 -0.0052199 -0.0005713  0.0052815  0.0185080 

Coefficients:
  Estimate Std. Error  t value Pr(|t|)
(Intercept) -1.3264070  0.0002227 -5957.00   2e-16 ***
s1  -1.3083551  0.0003149 -4154.91   2e-16 ***
s2  -0.0538174  0.0003149  -170.91   2e-16 ***
s3  -0.0346178  0.0003149  -109.94   2e-16 ***
s4  -0.0070427  0.0003149   -22.36   2e-16 ***
s5  -0.0069028  0.0003149   -21.92   2e-16 ***
s6   0.0155151  0.000314949.27   2e-16 ***
c1   0.3746976  0.0003149  1189.92   2e-16 ***
c2  -0.0149771  0.0003149   -47.56   2e-16 ***
c3   0.0272273  0.000314986.47   2e-16 ***
c4  -0.0282915  0.0003149   -89.84   2e-16 ***
c5  -0.0084926  0.0003149   -26.97   2e-16 ***
c6  -0.0093216  0.0003149   -29.60   2e-16 ***

Fundamental harmonic is there: a1 =  -1.3083551, b1 =   0.3746976 and, as
expected, higher order harmonics have much smaller coefficients. BUT
residuals are PERIODIC:

http://r.789695.n4.nabble.com/file/n4542389/residualsplot.jpg 

Since my error doesn't look like noise ( is strongly autocorrelated) did I
make any mistakes?

Thank you for reading down to this line...
I hope you have suggestions!

Thanks
Andrea


###
Here's a sample of data:

 data[1:10,] 
 time  forcevdef   hdef
1  0.  0.000   0.000  0.000
2  0.0025 -0.007  -0.229  0.000
3  0.0050 -0.166  -1.829  0.453
4  0.0075 -0.544  -7.086  2.038
5  0.0100 -0.957 -13.714  3.849
6  0.0125 -1.304 -19.886  5.660
7  0.0150 -1.550 -25.143  7.245
8  0.0175 -1.719 -29.714  8.377
9  0.0200 -1.838 -33.372  9.510
10 0.0225 -1.933 -36.343 10.642

--
View this message in context: 
http://r.789695.n4.nabble.com/Question-on-harmonic-Fourier-analysis-of-sinusoidal-time-series-tp4542389p4542389.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] Creating Better Table in R

2012-04-09 Thread bobo
Could anyone please direct me on how to make a nicer table in R? THANKS FOR
ALL THE HELP!

I would like to make a table with the following in it: estimate, t value,
significance, beta, standard errors, adjusted r squared, and residual
standard error (3 decimal points if possible, but I can do it by hand). Also
I'd love to include the source of my information on the bottom of the table
if possible.

file :  http://r.789695.n4.nabble.com/file/n4542349/datpat.csv datpat.csv 

this is the code I am using

mod.1-lm(Patents~FHouse, data=datpat)
summary(mod.1)
xtable(summary(mod.1))

*table I'm creating*
http://r.789695.n4.nabble.com/file/n4542349/chart.jpg 




*Detailed code*

*Summary Stats:*

Residuals:
Min  1Q  Median  3Q Max 
-1.7540 -0.8833 -0.5123  0.1183 11.9858 

Coefficients:
Estimate Std. Error t value Pr(|t|)
(Intercept)  2.617760.52866   4.952 1.34e-06 ***
FHouse  -0.187920.08489  -2.214   0.0277 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.77 on 254 degrees of freedom
Multiple R-squared: 0.01893,Adjusted R-squared: 0.01507 
F-statistic: 4.901 on 1 and 254 DF,  p-value: 0.02773 

*X Table Code:*
\begin{table}[ht]
\begin{center}
\begin{tabular}{r}
  \hline
  Estimate  Std. Error  t value  Pr($$$|$t$|$) \\ 
  \hline
(Intercept)  2.6178  0.5287  4.95  0. \\ 
  FHouse  -0.1879  0.0849  -2.21  0.0277 \\ 
   \hline
\end{tabular}
\end{center}
\end{table}

--
View this message in context: 
http://r.789695.n4.nabble.com/Creating-Better-Table-in-R-tp4542349p4542349.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] Shewhart Control Charts for Time Series data (qcc package)

2012-04-09 Thread Alexandr M
Hello everybody!

I am trying to apply qcc function to time series data.

There is an example from documentation:

###
library(qcc)
data(pistonrings)
attach(pistonrings)
diameter - qcc.groups(diameter, sample)
qcc(diameter[1:25,], type=xbar)
detach(pistonrings)
###

So, we have 5 iteams at each row - sample size =5

But what should I do if I have just a one dimensional time series data?

###
out.m - aggregate(pistonrings$diameter, list(sample = pistonrings$sample),
FUN=mean)
qcc(out.m$x, type=xbar) - Error: group sizes must be larger than one
###

Probably in this case samples are a consecutive subsets of the data but
they are not grouped as in an example.

I would like to apply qcc() like rollapply(width = sample_size, by =
step) in such a way that 'width' would stands for a samples taken from an
input data.

Thank you in advance !

--
Kind regards, Alexandr

[[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] i am not getting time series prediction results?

2012-04-09 Thread sagarnikam123
i have attached file,
http://r.789695.n4.nabble.com/file/n4542543/1BHP_A.txt 1BHP_A.txt 
i just want to validate predicting results, i give less 10 digits for
prediction  compare results it with previous input data(input file)

 i-read.table(file.choose())  #attached file taking
 i-i$V1
 length(i)
[1] 44
 j-i[1:34]
 pre-predict(ar(i),n.ahead=10)
 ts.plot(ts(j),pre$pred,col=c(1,5))
 pre$pred
Time Series:
Start = 45 
End = 54 
Frequency = 1 
 [1] -0.01136364 -0.01136364 -0.01136364 -0.01136364 -0.01136364
 [6] -0.01136364 -0.01136364 -0.01136364 -0.01136364 -0.01136364

prediction gives similar values,i.e. straight line.
why this so?

--
View this message in context: 
http://r.789695.n4.nabble.com/i-am-not-getting-time-series-prediction-results-tp4542543p4542543.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 Sweave users collaborate with Word users?

2012-04-09 Thread Richard M. Heiberger
You might want to consider SWord, which provides similar facilities for the
Word and R
user.  Word-oriented co-authors can modify the Word part of the document
without
impacting the R part of the document.

SWord is by Thomas Baier tho...@statconn.com, author of the statconnDCOM
interface
that is underneath RExcel.  See rcom.univie.ac.at for information and
download and to
sign up on the rcom email list.

Rich

On Sat, Apr 7, 2012 at 4:54 PM, Alexander Shenkin ashen...@ufl.edu wrote:

 Hello All,

 I'm getting my workflow switched over to Sweave, which is very cool.
 However, I collaborate with folks (as many of you must as well) who use
 Word to Track Changes amongst a group while crafting a paper.  In the
 simplest case, there will just be two people (one Sweave user and one
 Word user) editing a paper.

 I'm wondering, how do Sweave users go about this?  I could convert a
 sweave file to a .docx easily enough via an intermediary pdf, rtf, html
 or otherwise.  However, once the file has been marked up with changes,
 the challenge is to migrate those (accepted) changes back to the sweave
 document.  Perhaps the most straightforward way is to manually
 back-propagate changes, but I imagine that could be a painstaking process.

 Ideally, I imagine a tool that puts invisible tags in the word document
 when it is originally produced from Sweave, and is then able to
 propagate changes back to that sweave file after markup.  I'd be
 pleasantly surprised if such a tool existed.

 Perhaps there are other ways of making this work.  Any thoughts are
 kindly appreciated!

 Thanks,
 Allie

 __
 R-help@r-project.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.


[R] Help using R 2.14.2

2012-04-09 Thread Kinshuk Nayak
 Sir, 
 I am a student in biostat and bioinformatics.  I am interested to use R for my 
research work related to codon usage analysis.This software was used in a 
publication entitled  Online synonymous codon usage analyses with the ade4 and 
seqinR packages and a weblink  
http://pbil.univ-lyon1.fr/datasets/charif04/ was given for readers to use the 
program. Now  I would like to use the webpage for codon usage analysis for my 
selected genome.However, the program is not working for any other genome other 
than its default data set.I can filter the genome to keep only sequences more 
than 300bp in a text file in fasta formata and I would like to run the 
program/weblink for this file.

Could you please let me know, how to use the program for my dataset on the 
weblink / a edited R script for running the program on R 2.14.2 for this set of 
sequences.

Sincerely
Kinshuk Chandra Nayak
Institute of Life Sciences
Department of Biotechnology,Govt India
Bhubaneswar-751023,
India
[[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] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread David Winsemius


On Apr 8, 2012, at 11:48 PM, Guaramy wrote:


Sorry but i didn't understand.


The path to understanding is through study and application of the  
advice in the Posting Guide.


--
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] How do Sweave users collaborate with Word users?

2012-04-09 Thread Paul Bivand
If you're considering SWord, you should remember that the licence is
not the normal R licence and in commercial use will require a
commercial licence. While some academic disciplines use Word etc, the
issue may be more common outside academia.

For those of us where such requirements involve a procurement process,
the need to purchase something when other users (and the
administration) are happy with their Word/Excel solutions may be an
insurmountable barrier.

Good luck

Paul Bivand
Centre for Economic and Social Inclusion (a non-profit organisation)
London



On 9 April 2012 13:23, Richard M. Heiberger r...@temple.edu wrote:
 You might want to consider SWord, which provides similar facilities for the
 Word and R
 user.  Word-oriented co-authors can modify the Word part of the document
 without
 impacting the R part of the document.

 SWord is by Thomas Baier tho...@statconn.com, author of the statconnDCOM
 interface
 that is underneath RExcel.  See rcom.univie.ac.at for information and
 download and to
 sign up on the rcom email list.

 Rich

 On Sat, Apr 7, 2012 at 4:54 PM, Alexander Shenkin ashen...@ufl.edu wrote:

 Hello All,

 I'm getting my workflow switched over to Sweave, which is very cool.
 However, I collaborate with folks (as many of you must as well) who use
 Word to Track Changes amongst a group while crafting a paper.  In the
 simplest case, there will just be two people (one Sweave user and one
 Word user) editing a paper.

 I'm wondering, how do Sweave users go about this?  I could convert a
 sweave file to a .docx easily enough via an intermediary pdf, rtf, html
 or otherwise.  However, once the file has been marked up with changes,
 the challenge is to migrate those (accepted) changes back to the sweave
 document.  Perhaps the most straightforward way is to manually
 back-propagate changes, but I imagine that could be a painstaking process.

 Ideally, I imagine a tool that puts invisible tags in the word document
 when it is originally produced from Sweave, and is then able to
 propagate changes back to that sweave file after markup.  I'd be
 pleasantly surprised if such a tool existed.

 Perhaps there are other ways of making this work.  Any thoughts are
 kindly appreciated!

 Thanks,
 Allie

 __
 R-help@r-project.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.

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


Re: [R] Stepwise procedure with force.in command

2012-04-09 Thread Peter Ehlers

On 2012-04-09 00:44, Hien Nguyen wrote:

Dear R-helpers,

I am trying to do a stepwise procedure in which I want to force some
variables in the model. I have searched around and it seems that only
leaps package allows to force the variable in the stepwise procedure. I
use the leaps package and use the regsubsets(lm1, force.in = 1, data) to
force 1 variable in the model. However, the force.in command only allow
me to force 1 variable in. I want to force several variables in the
procedure. How could I force them in?


Well, I don't know where you've searched, but the help page for
step() or for stepAIC() in pkg MASS tells you that there's a
scope= argument which should do what you want.

Peter Ehlers



Thanks a lot,

Hien

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


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


[R] R problem nls

2012-04-09 Thread Karen Vandepoel
Hi,

I will try to explain what it is I need to do, how far I am in doing it yet
and where my problem is:

I have a lot of x,y values I need to fit a non linear function through.
Subsequently, I need to find the intersection point of this fitted curve
with y=1.01

The problem is I have a lot of values so I want to be able to do it all at
once.

I already imported my excel file of the points I have to use.

The things I already have are following (not all of the data is visible
because otherwise the file would be too long for this email, but i just
showed part of my data so you could better understand my problem):


__
  X1242   X1242.5 X1243   X1243.5 X1244   X1244.5
X1245   X1245.5 X1246   X1246.5 X1247   X1247.5 X1248
X1248.5 X1249 =  names of each data set ( corresponds to wavelengths,
ranges from 400 to 2500 with 0.5 steps = a lot of points!)
18.14 0.9860316 0.9860272 0.9860203 0.9860121 0.9860044 0.9859994 0.9859971
0.9859976 0.985 0.9860035 0.9860069 0.9860103 0.9860128 0.9860151
0.9860178
15.8  0.9857134 0.9857106 0.9857063 0.9857011 0.9856958 0.9856917 0.9856887
0.9856874 0.9856880 0.9856893 0.9856906 0.9856919 0.9856922 0.9856916
0.9856912
13.77 0.9930109 0.9930015 0.9929921 0.9929833 0.9929765 0.9929714 0.9929674
0.9929637 0.9929603 0.9929569 0.9929533 0.9929501 0.9929469 0.9929440
0.9929416
9.03  0.9875374 0.9875321 0.9875242 0.9875140 0.9875024 0.9874921 0.9874840
0.9874793 0.9874780 0.9874802 0.9874835 0.9874869 0.9874877 0.9874857
0.9874801
6.14  0.9900554 0.9900465 0.9900376 0.9900286 0.9900204 0.9900122 0.9900032
0.9899924 0.9899802 0.9899669 0.9899544 0.9899445 0.9899372 0.9899333
0.9899317
4.27  1.0050327 1.0050242 1.0050175 1.0050129 1.0050107 1.0050101 1.0050094
1.0050070 1.0050025 1.0049959 1.0049885 1.0049812 1.0049746 1.0049691
1.0049647
2.77  0.9892697 0.9892585 0.9892454 0.9892311 0.9892164 0.9892030 0.9891906
0.9891796 0.9891704 0.9891624 0.9891550 0.9891480 0.9891401 0.9891320
0.9891235
1.52  0.9979284 0.9979430 0.9979548 0.9979644 0.9979739 0.9979850 0.9979984
0.9980137 0.9980312 0.9980498 0.9980691 0.9980897 0.9981105 0.9981323
0.9981542
These are my x-values and y-values
X1249.5 X1250   X1250.5 X1251   X1251.5 X1252
X1252.5 X1253   X1253.5 X1254   X1254.5 X1255   X1255.5
X1256   X1256.5
18.14 0.9860214 0.9860261 0.9860320 0.9860377 0.9860425 0.9860456 0.9860462
0.9860449 0.9860433 0.9860422 0.9860417 0.9860428 0.9860444 0.9860456
0.9860456
15.8  0.9856911 0.9856918 0.9856934 0.9856958 0.9856984 0.9857014 0.9857040
0.9857069 0.9857099 0.9857132 0.9857153 0.9857156 0.9857132 0.9857080
0.9857016
13.77 0.9929406 0.9929405 0.9929425 0.9929445 0.9929464 0.9929476 0.9929476
0.9929455 0.9929431 0.9929409 0.9929377 0.9929356 0.9929331 0.9929304
0.9929279
9.03  0.9874721 0.9874641 0.9874578 0.9874542 0.9874529 0.9874541 0.9874556
0.9874563 0.9874565 0.9874551 0.9874527 0.9874498 0.9874461 0.9874415
0.9874371
6.14  0.9899318 0.9899319 0.9899304 0.9899263 0.9899192 0.9899095 0.9898986
0.9898873 0.9898777 0.9898703 0.9898641 0.9898591 0.9898546 0.9898495
0.9898439
4.27  1.0049605 1.0049564 1.0049528 1.0049495 1.0049466 1.0049435 1.0049404
1.0049357 1.0049298 1.0049230 1.0049155 1.0049093 1.0049049 1.0049019
1.0049002
2.77  0.9891158 0.9891100 0.9891058 0.9891036 0.9891013 0.9890986 0.9890943
0.9890873 0.9890789 0.9890691 0.9890583 0.9890485 0.9890395 0.9890323
0.9890266
1.52  0.9981760 0.9981970 0.9982176 0.9982372 0.9982565 0.9982753 0.9982935
0.9983102 0.9983259 0.9983408 0.9983533 0.9983647 0.9983749 0.9983841
0.9983929
  X1257   X1257.5 X1258   X1258.5 X1259   X1259.5
X1260   X1260.5 X1261   X1261.5 X1262   X1262.5 X1263
X1263.5 X1264
18.14 0.9860445 0.9860437 0.9860438 0.9860467 0.9860527 0.9860613 0.9860705
0.9860782 0.9860830 0.9860844 0.9860832 0.9860806 0.9860774 0.9860746
0.9860716
15.8  0.9856955 0.9856925 0.9856930 0.9856967 0.9857032 0.9857098 0.9857162
0.9857213 0.9857248 0.9857268 0.9857283 0.9857298 0.9857314 0.9857340
0.9857374
13.77 0.9929253 0.9929242 0.9929234 0.9929239 0.9929254 0.9929276 0.9929302
0.9929321 0.9929329 0.9929328 0.9929315 0.9929294 0.9929275 0.9929257
0.9929243
9.03  0.9874330 0.9874313 0.9874312 0.9874335 0.9874375 0.9874418 0.9874468
0.9874510 0.9874547 0.9874577 0.9874602 0.9874623 0.9874640 0.9874659
0.9874680
6.14  0.9898392 0.9898366 0.9898360 0.9898381 0.9898417 0.9898455 0.9898484
0.9898489 0.9898469 0.9898432 0.9898392 0.9898363 0.9898354 0.9898371
0.9898402
4.27  1.0048982 1.0048966 1.0048942 1.0048923 1.0048914 1.0048917 1.0048929
1.0048939 1.0048933 1.0048910 1.0048865 1.0048803 1.0048742 1.0048693
1.0048658
2.77  0.9890227 0.9890210 0.9890197 0.9890184 0.9890164 0.9890143 0.9890120
0.9890102 0.9890088 0.9890087 0.9890089 0.9890086 0.9890073 0.9890056
0.9890024
1.52  0.9984034 0.9984168 0.9984328 0.9984514 0.9984707 0.9984899 0.9985069
0.9985223 0.9985357 0.9985479 0.9985593 0.9985694 0.9985778 

Re: [R] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread Guaramy
I read it believe me . The reason that a i post this is because i am making a
thesis and a i am having this problem for over 2 weeks.
I can´t solve it and its causing me real problems.

Sorry for any inconvenience


Thanks  

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-integrate-int-fn-lower-0-upper-Inf-evaluation-of-function-gave-a-result-of-wrong-length-tp4541250p4542677.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] Creating Better Table in R

2012-04-09 Thread Carlos Ortega
Hi,

Yes, please check package tables.

Regards,
Carlos Ortega
www.qualityexcellence.es

2012/4/9 bobo bleza...@gmail.com

 Could anyone please direct me on how to make a nicer table in R? THANKS FOR
 ALL THE HELP!

 I would like to make a table with the following in it: estimate, t value,
 significance, beta, standard errors, adjusted r squared, and residual
 standard error (3 decimal points if possible, but I can do it by hand).
 Also
 I'd love to include the source of my information on the bottom of the table
 if possible.

 file :  http://r.789695.n4.nabble.com/file/n4542349/datpat.csv datpat.csv

 this is the code I am using

 mod.1-lm(Patents~FHouse, data=datpat)
 summary(mod.1)
 xtable(summary(mod.1))

 *table I'm creating*
 http://r.789695.n4.nabble.com/file/n4542349/chart.jpg




 *Detailed code*

 *Summary Stats:*

 Residuals:
Min  1Q  Median  3Q Max
 -1.7540 -0.8833 -0.5123  0.1183 11.9858

 Coefficients:
Estimate Std. Error t value Pr(|t|)
 (Intercept)  2.617760.52866   4.952 1.34e-06 ***
 FHouse  -0.187920.08489  -2.214   0.0277 *
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Residual standard error: 1.77 on 254 degrees of freedom
 Multiple R-squared: 0.01893,Adjusted R-squared: 0.01507
 F-statistic: 4.901 on 1 and 254 DF,  p-value: 0.02773

 *X Table Code:*
 \begin{table}[ht]
 \begin{center}
 \begin{tabular}{r}
  \hline
   Estimate  Std. Error  t value  Pr($$$|$t$|$) \\
  \hline
 (Intercept)  2.6178  0.5287  4.95  0. \\
  FHouse  -0.1879  0.0849  -2.21  0.0277 \\
   \hline
 \end{tabular}
 \end{center}
 \end{table}

 --
 View this message in context:
 http://r.789695.n4.nabble.com/Creating-Better-Table-in-R-tp4542349p4542349.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.




-- 
Saludos,
Carlos Ortega
www.qualityexcellence.es

[[alternative HTML version deleted]]

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


Re: [R] R problem nls

2012-04-09 Thread Jean V Adams
Your formula could be simplified to
y ~ 1 + a * exp(-b * x)
Solve this equation for x
x = ln[(y - 1)/a]/b
Use this equation to find the intersection point at a given value of y.
For example, when
y = 1.01
x = ln(0.01/a)/b

Jean


Karen Vandepoel wrote on 04/09/2012 07:33:19 AM:

 Hi,
 
 I will try to explain what it is I need to do, how far I am in doing it 
yet
 and where my problem is:
 
 I have a lot of x,y values I need to fit a non linear function through.
 Subsequently, I need to find the intersection point of this fitted curve
 with y=1.01
 
 The problem is I have a lot of values so I want to be able to do it all 
at
 once.
 
 I already imported my excel file of the points I have to use.
 
 The things I already have are following (not all of the data is visible
 because otherwise the file would be too long for this email, but i just
 showed part of my data so you could better understand my problem):
 
 
 __
   X1242   X1242.5 X1243   X1243.5 X1244   X1244.5
 X1245   X1245.5 X1246   X1246.5 X1247   X1247.5 X1248
 X1248.5 X1249 =  names of each data set ( corresponds to 
wavelengths,
 ranges from 400 to 2500 with 0.5 steps = a lot of points!)
 18.14 0.9860316 0.9860272 0.9860203 0.9860121 0.9860044 0.9859994 
0.9859971
 0.9859976 0.985 0.9860035 0.9860069 0.9860103 0.9860128 0.9860151
 0.9860178
 15.8  0.9857134 0.9857106 0.9857063 0.9857011 0.9856958 0.9856917 
0.9856887
 0.9856874 0.9856880 0.9856893 0.9856906 0.9856919 0.9856922 0.9856916
 0.9856912
 13.77 0.9930109 0.9930015 0.9929921 0.9929833 0.9929765 0.9929714 
0.9929674
 0.9929637 0.9929603 0.9929569 0.9929533 0.9929501 0.9929469 0.9929440
 0.9929416
 9.03  0.9875374 0.9875321 0.9875242 0.9875140 0.9875024 0.9874921 
0.9874840
 0.9874793 0.9874780 0.9874802 0.9874835 0.9874869 0.9874877 0.9874857
 0.9874801
 6.14  0.9900554 0.9900465 0.9900376 0.9900286 0.9900204 0.9900122 
0.9900032
 0.9899924 0.9899802 0.9899669 0.9899544 0.9899445 0.9899372 0.9899333
 0.9899317
 4.27  1.0050327 1.0050242 1.0050175 1.0050129 1.0050107 1.0050101 
1.0050094
 1.0050070 1.0050025 1.0049959 1.0049885 1.0049812 1.0049746 1.0049691
 1.0049647
 2.77  0.9892697 0.9892585 0.9892454 0.9892311 0.9892164 0.9892030 
0.9891906
 0.9891796 0.9891704 0.9891624 0.9891550 0.9891480 0.9891401 0.9891320
 0.9891235
 1.52  0.9979284 0.9979430 0.9979548 0.9979644 0.9979739 0.9979850 
0.9979984
 0.9980137 0.9980312 0.9980498 0.9980691 0.9980897 0.9981105 0.9981323
 0.9981542
 These are my x-values and y-values
 X1249.5 X1250   X1250.5 X1251   X1251.5 X1252
 X1252.5 X1253   X1253.5 X1254   X1254.5 X1255   X1255.5
 X1256   X1256.5
 18.14 0.9860214 0.9860261 0.9860320 0.9860377 0.9860425 0.9860456 
0.9860462
 0.9860449 0.9860433 0.9860422 0.9860417 0.9860428 0.9860444 0.9860456
 0.9860456
 15.8  0.9856911 0.9856918 0.9856934 0.9856958 0.9856984 0.9857014 
0.9857040
 0.9857069 0.9857099 0.9857132 0.9857153 0.9857156 0.9857132 0.9857080
 0.9857016
 13.77 0.9929406 0.9929405 0.9929425 0.9929445 0.9929464 0.9929476 
0.9929476
 0.9929455 0.9929431 0.9929409 0.9929377 0.9929356 0.9929331 0.9929304
 0.9929279
 9.03  0.9874721 0.9874641 0.9874578 0.9874542 0.9874529 0.9874541 
0.9874556
 0.9874563 0.9874565 0.9874551 0.9874527 0.9874498 0.9874461 0.9874415
 0.9874371
 6.14  0.9899318 0.9899319 0.9899304 0.9899263 0.9899192 0.9899095 
0.9898986
 0.9898873 0.9898777 0.9898703 0.9898641 0.9898591 0.9898546 0.9898495
 0.9898439
 4.27  1.0049605 1.0049564 1.0049528 1.0049495 1.0049466 1.0049435 
1.0049404
 1.0049357 1.0049298 1.0049230 1.0049155 1.0049093 1.0049049 1.0049019
 1.0049002
 2.77  0.9891158 0.9891100 0.9891058 0.9891036 0.9891013 0.9890986 
0.9890943
 0.9890873 0.9890789 0.9890691 0.9890583 0.9890485 0.9890395 0.9890323
 0.9890266
 1.52  0.9981760 0.9981970 0.9982176 0.9982372 0.9982565 0.9982753 
0.9982935
 0.9983102 0.9983259 0.9983408 0.9983533 0.9983647 0.9983749 0.9983841
 0.9983929
   X1257   X1257.5 X1258   X1258.5 X1259   X1259.5
 X1260   X1260.5 X1261   X1261.5 X1262   X1262.5 X1263
 X1263.5 X1264
 18.14 0.9860445 0.9860437 0.9860438 0.9860467 0.9860527 0.9860613 
0.9860705
 0.9860782 0.9860830 0.9860844 0.9860832 0.9860806 0.9860774 0.9860746
 0.9860716
 15.8  0.9856955 0.9856925 0.9856930 0.9856967 0.9857032 0.9857098 
0.9857162
 0.9857213 0.9857248 0.9857268 0.9857283 0.9857298 0.9857314 0.9857340
 0.9857374
 13.77 0.9929253 0.9929242 0.9929234 0.9929239 0.9929254 0.9929276 
0.9929302
 0.9929321 0.9929329 0.9929328 0.9929315 0.9929294 0.9929275 0.9929257
 0.9929243
 9.03  0.9874330 0.9874313 0.9874312 0.9874335 0.9874375 0.9874418 
0.9874468
 0.9874510 0.9874547 0.9874577 0.9874602 0.9874623 0.9874640 0.9874659
 0.9874680
 6.14  0.9898392 0.9898366 0.9898360 0.9898381 0.9898417 0.9898455 
0.9898484
 0.9898489 0.9898469 0.9898432 0.9898392 0.9898363 0.9898354 0.9898371
 0.9898402
 4.27  1.0048982 

Re: [R] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread Berend Hasselman

On 09-04-2012, at 14:54, Guaramy wrote:

 I read it believe me . The reason that a i post this is because i am making a
 thesis and a i am having this problem for over 2 weeks.
 I can´t solve it and its causing me real problems.


You didn't give us a reproducible example.
I'm assuming your original function is

GNL.pdf.fn - function(x,mu,sigma,alpha,beta,rho) # missing in your original 
post
{
y  -  x-rho*mu

cf.fn  -  function(s){
cplex = complex(1,0,1)
temp1 = alpha*beta*exp(-sigma*s^2/2)
temp2 = (alpha-cplex*s)*(beta+cplex*s)
out = (temp1/temp2)^rho
out
}

temp.fn  -  function(s){ Mod(cf.fn(s))*cos(Arg(cf.fn(s))-s*y) }

int.fn  -  function(t){sapply(t,FUN=temp.fn)}
te  -  
integrate(int.fn,lower=0,upper=Inf,rel.tol=1e-10,subdivisions=100)
temp3  -  ifelse(te$message == OK,te$value/pi,NA)
temp3
}

Question: why are you using ifelse when te$message is a scalar?

if( te$message == OK ) te$value/pi else  NA

would do what you want.

David meant that you should have a look at Vectorize. And then do some 
experimenting.
So I tried this

GNL.pdf.fn(1, .2, .3, 1, 1, .6)

GNL.vec - Vectorize(GNL.pdf.fn, x)

GNL.vec(c(.9,1,1.1), .2, .3, 1, 1, .6)

And now it's your turn for seq(-4,4,0.1).

Berend

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


Re: [R] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread peter dalgaard

On Apr 9, 2012, at 14:54 , Guaramy wrote:

 I read it believe me . The reason that a i post this is because i am making a
 thesis and a i am having this problem for over 2 weeks.
 I can´t solve it and its causing me real problems.

Well, if you have a function that is not vectorized, i.e. it works for a single 
argument but not for a vector of arguments, and you want it to be vectorized, 
the tool is Vectorize(), the help page of which David pointed you to.  (This 
essentially works around the problem by creating a new function which calls the 
function one vector element at the time.)

The direct cause of your problem seems to be that if x is a vector, so is y, 
and your temp.fn returns a vector, so int.fn returns a matrix, and integrate() 
is unhappy with that.

-- 
Peter Dalgaard, Professor,
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.


[R] Specifying the ordering of a vector

2012-04-09 Thread crmnaw
Hi, 

I'm trying to create a vector (or matrix row) with a specific ordering. For
example, I have the following vector:

x-c(0.1,0.2,0.3,0.4,0.5,0.6) 

that has order

order(x)
[1] 1 2 3 4 5 6

I want another vector that has the same values as x, but with a different
ordering. For example, I want y to have values 0.1, 0.2, etc. but in the
order 1-2-5-6-3-4. The answer would be

y
[1] 0.1 0.2 0.5 0.6 0.3 0.4

Any help would be greatly appreciated. Thanks. 

--
View this message in context: 
http://r.789695.n4.nabble.com/Specifying-the-ordering-of-a-vector-tp4542766p4542766.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] Comparing 2 means. pls help

2012-04-09 Thread R. Michael Weylandt
You should look up what a t-test is.

Michael

On Mon, Apr 9, 2012 at 2:58 AM, ali_protocol
mohammadianalimohammad...@gmail.com wrote:
 I am interested in the difference of 2 data:

 mat1= c(2.2, 2.3, 2.2,2.5)
 mat2= c(2.6, 2.8, 2.7,2.4)

 mat= mat2-mat1

 I perform an action on both mat1 and mat2, and I get
 mat1prime and mat2prime:

 mat1prime= c(2.5, 2.5, 2.3,2.5)
 mat2prime= c(2.6, 2.8, 2.7,2.6)

 matprime= mat2prime-mat1prime


 I want to check wether mat and matprime have the same mean.
 How should I do this in R?

 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Comparing-2-means-pls-help-tp4542237p4542237.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] Avoid loop with the integrate function

2012-04-09 Thread Navin Goyal
Hi Berend,
I now understand what you were explaining. It works great.
Thanks a lot for your time and help. I really appreciate it.

Navin Goyal



On Mon, Apr 9, 2012 at 1:26 AM, Berend Hasselman b...@xs4all.nl wrote:


 On 09-04-2012, at 01:26, Navin Goyal wrote:

  Hi,
  I am not sure I follow the scalar part of the suggestion. Could you
 please elaborate it a bit more?  Each row has a different value of score
 for each subject at each time point.

 That last remark has nothing to do with it.
 The point is that in your function func1 the returned value doesn't depend
 on the argument t.
 So func1 is returning a vector of length==length(t) of identical values
 except for the first entry which is zero.

 Do this after the loop for comb3


 # No integral since expression in func1 doesn't depend on t
 comb4-comb2
 comb4$cmhz=0
 comb4-comb4[order(comb4$ID, comb4$TIME), ]
 for (q in 1:length(comb4$ID))
 {
comb4$cmhz[q]-comb4$TIME[q]*func1(comb4$TIME[q],
cov1=0.001,beta1=0.02,
change=comb4$score[q], other=comb4$frac[q])
 }
 head(comb4)

 all.equal(comb3$cmhz, comb4$cmhz)

 You don't need integrate to get what you seem to be wanting.

 Berend

  Also I simplified the example but there are other terms that  change
 over each row for each subject.
  Does this mean that loops is the only way to go ?
 
  Thanks again for your help and time.
 
  Navin Goyal
 
 
  On Sun, Apr 8, 2012 at 4:03 AM, Berend Hasselman b...@xs4all.nl wrote:
 
  On 08-04-2012, at 08:28, Navin Goyal wrote:
 
   Dear R users,
   I am running a loop with the integrate function. I have pasted the code
   below. I am integrating a function from time=0 to the time value in
 every
   row.
   I have to perform this integration over thousands of rows with
 different
   parameters in each row. Could someone please suggest if there is an
   efficient/faster/easier way to do this by avoiding the loops ?
  
   Thank you so much for your help and time.
   --
   Navin Goyal
  
   #
   dose-10
   time-0:5
   id-1:5
   data1-expand.grid(id,time,dose)
   names(data1)-c(ID,TIME, DOSE)
   data1-data1[order(data1$ID,data1$TIME),]
  
   ed-data1[!duplicated(data1$ID) , c(ID,DOSE)]
   set.seed(5324123)
  
   for (k in 1:length(ed$ID))
   {
   ed$base[k]-100*exp(rnorm(1,0,0.05))
   ed$drop[k]-0.2*exp(rnorm(1,0,0.01))
   ed$frac[k]-0.5*exp(rnorm(1,0,0.1))
   }
 
  Why not
 
  ed$base - 100*exp(rnorm(length(ed$ID), 0, 0.05))
 
  etc.
 
   comb1-merge(data1[, c(ID,TIME)], ed)
   comb2-comb1
   comb2$score-comb2$base*exp(-comb2$drop*comb2$TIME)
  
   func1-function(t,cov1,beta1, change,other)
   {
   ifelse(t==0,cov1, cov1*exp(beta1*change+other))
   }
 
  AFAICS, cov1, beta1, change and other are scalars.
  So when an item of t is == 0 then the function value is cov1 and in all
 other cases it is the same scalar and does not depend on t. So func1 is a
 step function. You could calculate the integral directly.
 
  Berend
 
 
 
 
  --
  Navin Goyal
 




-- 
Navin Goyal

[[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] Specifying the ordering of a vector

2012-04-09 Thread Duncan Murdoch

On 09/04/2012 9:38 AM, crmnaw wrote:

Hi,

I'm trying to create a vector (or matrix row) with a specific ordering. For
example, I have the following vector:

x-c(0.1,0.2,0.3,0.4,0.5,0.6)

that has order

order(x)
[1] 1 2 3 4 5 6

I want another vector that has the same values as x, but with a different
ordering. For example, I want y to have values 0.1, 0.2, etc. but in the
order 1-2-5-6-3-4. The answer would be

y
[1] 0.1 0.2 0.5 0.6 0.3 0.4

Any help would be greatly appreciated. Thanks.


x[c(1,2,5,6,3,4)]

will do it.

Duncan Murdoch

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


Re: [R] Specifying the ordering of a vector

2012-04-09 Thread crmnaw
Hi,

Thanks for your reply. In the example I gave in the original post, your code
works. But for others, it doesn't and I'm not sure why it works for some
cases and not for others. For example:

x-c(0.04,0.07,0.20,0.35,0.55,0.70)
order(x)
[1] 1 2 3 4 5 6 

Let's say I want y to be in the order 1-2-5-3-6-4. So I do:

y-x[c(1,2,5,3,6,4)]
y
[1] 0.04 0.07 0.55 0.20 0.70 0.35
order(y)
[1] 1 2 4 6 3 5

So, y isn't in the order I want it to be in. Any help that can be provided
would be greatly appreciated. 


--
View this message in context: 
http://r.789695.n4.nabble.com/Specifying-the-ordering-of-a-vector-tp4542766p4542815.html
Sent from the R help mailing list archive at Nabble.com.

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


Re: [R] R problem nls

2012-04-09 Thread Carlos Ortega
Hi,

And regarding how to extend the nls algorithm to a larger dataset it is a
question of indicating in the nls() functions which data.frame to use.
In you example you were using a small set of x's ad y's, so if you want to
use a large set put it in a new data.frame and pass it to nls().

And if you want to repeat that calculation for many different sets of x,y
(each one of them corresponding to a different wavelength) you could do
that in a loop.

Regards,
Carlos Ortega
www.qualityexcellence.es

2012/4/9 Karen Vandepoel karen.vandep...@gmail.com

 Hi,

 I will try to explain what it is I need to do, how far I am in doing it yet
 and where my problem is:

 I have a lot of x,y values I need to fit a non linear function through.
 Subsequently, I need to find the intersection point of this fitted curve
 with y=1.01

 The problem is I have a lot of values so I want to be able to do it all at
 once.

 I already imported my excel file of the points I have to use.

 The things I already have are following (not all of the data is visible
 because otherwise the file would be too long for this email, but i just
 showed part of my data so you could better understand my problem):


 __
  X1242   X1242.5 X1243   X1243.5 X1244   X1244.5
 X1245   X1245.5 X1246   X1246.5 X1247   X1247.5 X1248
 X1248.5 X1249 =  names of each data set ( corresponds to wavelengths,
 ranges from 400 to 2500 with 0.5 steps = a lot of points!)
 18.14 0.9860316 0.9860272 0.9860203 0.9860121 0.9860044 0.9859994 0.9859971
 0.9859976 0.985 0.9860035 0.9860069 0.9860103 0.9860128 0.9860151
 0.9860178
 15.8  0.9857134 0.9857106 0.9857063 0.9857011 0.9856958 0.9856917 0.9856887
 0.9856874 0.9856880 0.9856893 0.9856906 0.9856919 0.9856922 0.9856916
 0.9856912
 13.77 0.9930109 0.9930015 0.9929921 0.9929833 0.9929765 0.9929714 0.9929674
 0.9929637 0.9929603 0.9929569 0.9929533 0.9929501 0.9929469 0.9929440
 0.9929416
 9.03  0.9875374 0.9875321 0.9875242 0.9875140 0.9875024 0.9874921 0.9874840
 0.9874793 0.9874780 0.9874802 0.9874835 0.9874869 0.9874877 0.9874857
 0.9874801
 6.14  0.9900554 0.9900465 0.9900376 0.9900286 0.9900204 0.9900122 0.9900032
 0.9899924 0.9899802 0.9899669 0.9899544 0.9899445 0.9899372 0.9899333
 0.9899317
 4.27  1.0050327 1.0050242 1.0050175 1.0050129 1.0050107 1.0050101 1.0050094
 1.0050070 1.0050025 1.0049959 1.0049885 1.0049812 1.0049746 1.0049691
 1.0049647
 2.77  0.9892697 0.9892585 0.9892454 0.9892311 0.9892164 0.9892030 0.9891906
 0.9891796 0.9891704 0.9891624 0.9891550 0.9891480 0.9891401 0.9891320
 0.9891235
 1.52  0.9979284 0.9979430 0.9979548 0.9979644 0.9979739 0.9979850 0.9979984
 0.9980137 0.9980312 0.9980498 0.9980691 0.9980897 0.9981105 0.9981323
 0.9981542
 These are my x-values and y-values
X1249.5 X1250   X1250.5 X1251   X1251.5 X1252
 X1252.5 X1253   X1253.5 X1254   X1254.5 X1255   X1255.5
 X1256   X1256.5
 18.14 0.9860214 0.9860261 0.9860320 0.9860377 0.9860425 0.9860456 0.9860462
 0.9860449 0.9860433 0.9860422 0.9860417 0.9860428 0.9860444 0.9860456
 0.9860456
 15.8  0.9856911 0.9856918 0.9856934 0.9856958 0.9856984 0.9857014 0.9857040
 0.9857069 0.9857099 0.9857132 0.9857153 0.9857156 0.9857132 0.9857080
 0.9857016
 13.77 0.9929406 0.9929405 0.9929425 0.9929445 0.9929464 0.9929476 0.9929476
 0.9929455 0.9929431 0.9929409 0.9929377 0.9929356 0.9929331 0.9929304
 0.9929279
 9.03  0.9874721 0.9874641 0.9874578 0.9874542 0.9874529 0.9874541 0.9874556
 0.9874563 0.9874565 0.9874551 0.9874527 0.9874498 0.9874461 0.9874415
 0.9874371
 6.14  0.9899318 0.9899319 0.9899304 0.9899263 0.9899192 0.9899095 0.9898986
 0.9898873 0.9898777 0.9898703 0.9898641 0.9898591 0.9898546 0.9898495
 0.9898439
 4.27  1.0049605 1.0049564 1.0049528 1.0049495 1.0049466 1.0049435 1.0049404
 1.0049357 1.0049298 1.0049230 1.0049155 1.0049093 1.0049049 1.0049019
 1.0049002
 2.77  0.9891158 0.9891100 0.9891058 0.9891036 0.9891013 0.9890986 0.9890943
 0.9890873 0.9890789 0.9890691 0.9890583 0.9890485 0.9890395 0.9890323
 0.9890266
 1.52  0.9981760 0.9981970 0.9982176 0.9982372 0.9982565 0.9982753 0.9982935
 0.9983102 0.9983259 0.9983408 0.9983533 0.9983647 0.9983749 0.9983841
 0.9983929
  X1257   X1257.5 X1258   X1258.5 X1259   X1259.5
 X1260   X1260.5 X1261   X1261.5 X1262   X1262.5 X1263
 X1263.5 X1264
 18.14 0.9860445 0.9860437 0.9860438 0.9860467 0.9860527 0.9860613 0.9860705
 0.9860782 0.9860830 0.9860844 0.9860832 0.9860806 0.9860774 0.9860746
 0.9860716
 15.8  0.9856955 0.9856925 0.9856930 0.9856967 0.9857032 0.9857098 0.9857162
 0.9857213 0.9857248 0.9857268 0.9857283 0.9857298 0.9857314 0.9857340
 0.9857374
 13.77 0.9929253 0.9929242 0.9929234 0.9929239 0.9929254 0.9929276 0.9929302
 0.9929321 0.9929329 0.9929328 0.9929315 0.9929294 0.9929275 0.9929257
 0.9929243
 9.03  0.9874330 0.9874313 0.9874312 0.9874335 0.9874375 0.9874418 0.9874468
 0.9874510 0.9874547 0.9874577 0.9874602 0.9874623 

[R] Calculating the overlapping area of ellipsoides

2012-04-09 Thread Steve_Friedman

I'm generating several scatter3d visualizations of vegetation community
distribution against a number parameters using the scatter3d function car
package,

The visualizations are fantastic and they provide great support for the
analysis. As far as I can determine, scatter3d does not save the output of
the routine as an object so I can not


x -  scatter3d(MADep2010  ~ Day7Max2010 + S7DryFreq2010  |
as.factor(VegClass), data = Veg ,
axis.scale = TRUE, sphere.size=1.5, surface=FALSE, parallel=FALSE,
ellipsoid=TRUE, revolutions=1,
   surface.col =c(red, green, blue, gold, magenta))

str(x)
NULL

I would like to extend the work by calculating the percent of overlapping
area among the ellipsoides on a community by community basis.

Can anyone point me to a package or suggest an approach to accomplishing
this.

Thanks
Steve


Steve Friedman Ph. D.
Ecologist  / Spatial Statistical Analyst
Everglades and Dry Tortugas National Park
950 N Krome Ave (3rd Floor)
Homestead, Florida 33034

steve_fried...@nps.gov
Office (305) 224 - 4282
Fax (305) 224 - 4147

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


Re: [R] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread Guaramy
Tanks a lot to all of you that take the time to help me.
This is a really useful and helpful forum  and i will try to help as much as
i can

Best Regards Jorge 

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-integrate-int-fn-lower-0-upper-Inf-evaluation-of-function-gave-a-result-of-wrong-length-tp4541250p4542842.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] Error in integrate(int.fn, lower = 0, upper = Inf) : evaluation of function gave a result of wrong length

2012-04-09 Thread Uwe Ligges



On 09.04.2012 16:08, Guaramy wrote:

Tanks a lot to all of you that take the time to help me.
This is a really useful and helpful forum  and i will try to help as much as
i can


Thanks for offering help to the R-help *mailing list*.
Please try to use the mailing list interface (rather than Nabble) and 
quote the context. Also, please send your message not only to the list, 
but also to those who need helped (or in this case, who helped you).


Best,
Uwe Ligges


Best Regards Jorge

--
View this message in context: 
http://r.789695.n4.nabble.com/Error-in-integrate-int-fn-lower-0-upper-Inf-evaluation-of-function-gave-a-result-of-wrong-length-tp4541250p4542842.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] Specifying the ordering of a vector

2012-04-09 Thread Duncan Murdoch

On 09/04/2012 9:58 AM, crmnaw wrote:

Hi,

Thanks for your reply. In the example I gave in the original post, your code
works. But for others, it doesn't and I'm not sure why it works for some
cases and not for others. For example:

x-c(0.04,0.07,0.20,0.35,0.55,0.70)
order(x)
[1] 1 2 3 4 5 6

Let's say I want y to be in the order 1-2-5-3-6-4. So I do:

y-x[c(1,2,5,3,6,4)]
y
[1] 0.04 0.07 0.55 0.20 0.70 0.35
order(y)
[1] 1 2 4 6 3 5

So, y isn't in the order I want it to be in. Any help that can be provided
would be greatly appreciated.


Okay, I didn't understand what you meant.

You should think of the result of order(y) as a permutation to apply to 
y so that it is in sorted order.  That is,


y[order(y)]

will always be in sorted order.  So what you want is a permutation of x 
that will be put into sorted order by y[c(1,2,5,3,6,4)], i.e. you want 
the inverse of the permutation c(1,2,5,3,6,4).  Turns out that if you 
apply order to a permutation you get its inverse.  So what you want is


y - x[order(c(1,2,5,3,6,4))]

Duncan Murdoch

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


Re: [R] How do Sweave users collaborate with Word users?

2012-04-09 Thread Alexander Shenkin
Thanks Rich,

While it doesn't tickle me the way sweave/knitr does, SWord sounds more
or less like the thing I'm looking for.  However, poking around
http://rcom.univie.ac.at/ and http://www.statconn.com doesn't reveal any
download links.  It seems as though SWord has been pulled from their
lineup?  Not sure - I've shot them an email to inquire.

Thanks,
Allie

On 4/9/2012 7:23 AM, Richard M. Heiberger wrote:
 You might want to consider SWord, which provides similar facilities
 for the Word and R
 user.  Word-oriented co-authors can modify the Word part of the
 document without
 impacting the R part of the document.
  
 SWord is by Thomas Baier tho...@statconn.com
 mailto:tho...@statconn.com, author of the statconnDCOM interface
 that is underneath RExcel.  See rcom.univie.ac.at
 http://rcom.univie.ac.at for information and download and to
 sign up on the rcom email list.
  
 Rich

 On Sat, Apr 7, 2012 at 4:54 PM, Alexander Shenkin ashen...@ufl.edu
 mailto:ashen...@ufl.edu wrote:

 Hello All,

 I'm getting my workflow switched over to Sweave, which is very cool.
 However, I collaborate with folks (as many of you must as well)
 who use
 Word to Track Changes amongst a group while crafting a paper.  In the
 simplest case, there will just be two people (one Sweave user and one
 Word user) editing a paper.

 I'm wondering, how do Sweave users go about this?  I could convert a
 sweave file to a .docx easily enough via an intermediary pdf, rtf,
 html
 or otherwise.  However, once the file has been marked up with changes,
 the challenge is to migrate those (accepted) changes back to the
 sweave
 document.  Perhaps the most straightforward way is to manually
 back-propagate changes, but I imagine that could be a painstaking
 process.

 Ideally, I imagine a tool that puts invisible tags in the word
 document
 when it is originally produced from Sweave, and is then able to
 propagate changes back to that sweave file after markup.  I'd be
 pleasantly surprised if such a tool existed.

 Perhaps there are other ways of making this work.  Any thoughts are
 kindly appreciated!

 Thanks,
 Allie

 __
 R-help@r-project.org mailto:R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 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] How do Sweave users collaborate with Word users?

2012-04-09 Thread Alexander Shenkin
Thanks for the heads up, Paul.  That's good to know.  I happen to be in
academics, as are my collaborators, but I could imagine running into
problems down the road.  The good thing seems to be that it's just the
users who want to interact with R who need the software.  If
collaborators are just touching the text, or are inserting graphics of
their own, then sounds like things are still good.  But, that's just my
impression from the brief description I've seen online.

Thanks,
Allie

On 4/9/2012 7:42 AM, Paul Bivand wrote:
 If you're considering SWord, you should remember that the licence is
 not the normal R licence and in commercial use will require a
 commercial licence. While some academic disciplines use Word etc, the
 issue may be more common outside academia.

 For those of us where such requirements involve a procurement process,
 the need to purchase something when other users (and the
 administration) are happy with their Word/Excel solutions may be an
 insurmountable barrier.

 Good luck

 Paul Bivand
 Centre for Economic and Social Inclusion (a non-profit organisation)
 London



 On 9 April 2012 13:23, Richard M. Heiberger r...@temple.edu wrote:
 You might want to consider SWord, which provides similar facilities for the
 Word and R
 user.  Word-oriented co-authors can modify the Word part of the document
 without
 impacting the R part of the document.

 SWord is by Thomas Baier tho...@statconn.com, author of the statconnDCOM
 interface
 that is underneath RExcel.  See rcom.univie.ac.at for information and
 download and to
 sign up on the rcom email list.

 Rich

 On Sat, Apr 7, 2012 at 4:54 PM, Alexander Shenkin ashen...@ufl.edu wrote:

 Hello All,

 I'm getting my workflow switched over to Sweave, which is very cool.
 However, I collaborate with folks (as many of you must as well) who use
 Word to Track Changes amongst a group while crafting a paper.  In the
 simplest case, there will just be two people (one Sweave user and one
 Word user) editing a paper.

 I'm wondering, how do Sweave users go about this?  I could convert a
 sweave file to a .docx easily enough via an intermediary pdf, rtf, html
 or otherwise.  However, once the file has been marked up with changes,
 the challenge is to migrate those (accepted) changes back to the sweave
 document.  Perhaps the most straightforward way is to manually
 back-propagate changes, but I imagine that could be a painstaking process.

 Ideally, I imagine a tool that puts invisible tags in the word document
 when it is originally produced from Sweave, and is then able to
 propagate changes back to that sweave file after markup.  I'd be
 pleasantly surprised if such a tool existed.

 Perhaps there are other ways of making this work.  Any thoughts are
 kindly appreciated!

 Thanks,
 Allie

 __
 R-help@r-project.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.
 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/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] i am not getting time series prediction results?

2012-04-09 Thread Uwe Ligges



On 09.04.2012 13:27, sagarnikam123 wrote:

i have attached file,
http://r.789695.n4.nabble.com/file/n4542543/1BHP_A.txt 1BHP_A.txt
i just want to validate predicting results, i give less 10 digits for
prediction  compare results it with previous input data(input file)


i-read.table(file.choose())  #attached file taking
i-i$V1
length(i)

[1] 44

j-i[1:34]
pre-predict(ar(i),n.ahead=10)
ts.plot(ts(j),pre$pred,col=c(1,5))
pre$pred

Time Series:
Start = 45
End = 54
Frequency = 1
  [1] -0.01136364 -0.01136364 -0.01136364 -0.01136364 -0.01136364
  [6] -0.01136364 -0.01136364 -0.01136364 -0.01136364 -0.01136364

prediction gives similar values,i.e. straight line.
why this so?


 ar(i)

gives:

Call:
ar(x = i)


Order selected 0  sigma^2 estimated as  0.0

hence an AR 0 process was fitted, i.e. no ar component, just th mean is 
used for prediction (and is the only estimated parameter).


 mean(i)
 -0.01136364

Uwe Ligges








--
View this message in context: 
http://r.789695.n4.nabble.com/i-am-not-getting-time-series-prediction-results-tp4542543p4542543.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] Difference between spec.pgram spec.ar

2012-04-09 Thread Uwe Ligges



On 08.04.2012 20:39, Bazman76 wrote:

Hi there,

Can someone explain what the difference between spec.pgram and spec.ar is?

I understand that they attempt to do the same thing one using an AR
estimation of the underlying series to estimate teh sensity the other using
the FFT. However when applied to teh same data set they seem to be giving
quite different results?

http://r.789695.n4.nabble.com/file/n4541358/R_example.jpg


Clearly the spec.ar() result seems to be smoothed but the results are also
generally very different only really sharing the peak as the frequencies go
towards zero.

Can someone please explain why these two functions produce such different
results on the same data set?


Because they really measure different things? Why do you expect to get 
the same output in time as well as in frequency domain?


Uwe Ligges




code shown below:

  library(waveslim)

vols=read.csv(file=C:/Users/ocuk/My Documents/Abs Vol.csv, header=TRUE,
sep=,)
x-ts(vols[,1])
#x

## LA(8)
vol.la8- mra(x, la8, 4, modwt)
names(vol.la8)- c(d1, d2, d3, d4, s4)
## plot multiresolution analysis of IBM data
#par(mfcol=c(6,1), pty=m, mar=c(5-2,4,4-2,2))
par(mfcol=c(3,1), pty=m, mar=c(5-2,4,4-2,2))
plot.ts(x, axes=F, ylab=, main=(abs rtns))
#for(i in 1:5)
#   plot.ts(vol.la8[[i]], axes=F, ylab=names(vol.la8)[i])
#axis(side=1, at=seq(0,1600,by=100),
#   
labels=c(0,,200,,400,,600,,800,,1000,,1200,,1400,,1600))

spectrum(vol.la8[[1]])
specx- spec.ar(x)

--
View this message in context: 
http://r.789695.n4.nabble.com/Difference-between-spec-pgram-spec-ar-tp4541358p4541358.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] Panel.abline would not show beyond a certain slope value

2012-04-09 Thread Uwe Ligges



On 09.04.2012 06:31, Aziz, Muhammad Fayez wrote:


Hi,

Please find the code and data following. Problem appears in lines 37 - 42 of 
the code. I am using R.2.13.0 on WinXP.


Works for me with a recent version of R (i.e. 2.15.0) / grid and lattice 
(i.e. 0.20-6).

Please always try with recent versions of R and packages before asking.

Best,
Uwe Ligges




Regards,
Fayez


PowerLawGraphsAblineQn.R‎

library('reshape')
library(lattice)
library(igraph) # to use power.law.fit function
library(latticeExtra) # to use panel.lmlineq in loglog xyplot

File2Open = C:\\Documents and Settings\\All 
Users\\Documents\\Desktop\\Fayez\\RPractice\\PowerLawGraphsAblineQnData.txt
DataTable = read.table(File2Open, header = TRUE, sep = \t)
md- melt(DataTable, id.vars = c('Range', 'TheValue', 'TheNo')) # removed 
'DegType' column 030312
md$variable- as.numeric(substr(md$variable, 9, 
nchar(as.character(md$variable


mypanel4loglog-
 function(x, # x is the variable column in melted data, equals the Domain 
Nos
 y, # y is the value column in melted data, the degrees
 ... # Rest of the arguments
 )
{

 kfreq- table(y); # compute frquency hash table of y, the values
 k- 1:max(y)
 for (i in k) {
 ichar- as.character(i) # convert to match the names(freq), the 
character-based hash key of freq, which is the value
 if (!(ichar %in% names(kfreq)))
 kfreq[ichar]- 0
 }
 sortedkeys- as.character(k)
 kfreq-  kfreq[sortedkeys]
 pk- kfreq / length(y)
 logk- log(k)
 logpk- log(pk)

 logpk[logpk == -Inf] =  # remove the -Inf or log(p(k)) = 0 values for lm function, 
NULL is 0-length so use  instead that has length of one null character
 logpk- as.numeric(logpk) #  converts all values to character, lm needs 
numeric
 print(rbind(logk, logpk))#write.table(rbind(k, kfreq, pk, logk, logpk), paste(FilePath, 
\\data, sep=), sep = \t, append = TRUE)

 panel.xyplot(col=blue, logk, logpk, type = c('p', 'r'))

 lm.r = lm(logpk ~ logk)
 panel.abline(coef=c(-4.847634, -1.037480)) # ---  This gets drawn
 panel.abline(coef=c(-4.847634, -1.037481)) # ---  This doesn't get drawn
 print(coef(lm.r)) # -4.847634   -1.349699 ;  -3.377894   -1.498693

} # end mypanel4loglog

myprepanel4loglog-
 function(x, # x is the variable column in melted data, equals the 
Network's ages
 y, # y is the value column in melted data, the degrees
 ... # Rest of the arguments
 )
{
 FreqTable- as.data.frame(table(y))
 FreqsVector- sort(FreqTable$Freq)
 Min- FreqsVector[1] # first element - the lowest value frequency
#print(c(Max2, length(y)))
 list(ylim = c(log(Min / length(y)), 0), xlim = c(0, log(max(y # 
log(p(k)) is always -ve as p(k) is decimal, so max(log(p(k)) is 0
} # end myprepanel4loglog

print(xyplot(value ~ variable | Range,
 data = md,
 xlab = log(k); Panel = Range,
 ylab = log(p(k)),
 main = log(k) vs. log(p(k)),
 groups = Range,
 pch = 20, # dots instead of circles
 panel = mypanel4loglog,
 prepanel = myprepanel4loglog, # to set the scale of k and pk
 scales = list(x = (relation = free), y = (relation = free)), # 
necessary to make x-axis in each panel adjustable according to k
))




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


Re: [R] Specifying the ordering of a vector

2012-04-09 Thread crmnaw
Thank you very much! This worked. 

--
View this message in context: 
http://r.789695.n4.nabble.com/Specifying-the-ordering-of-a-vector-tp4542766p4542907.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] Creating a dummy variable

2012-04-09 Thread Uwe Ligges



On 09.04.2012 07:16, arunkumar wrote:

HI

I want to create a dummy variable for a dataset which has various dates Eg
I've one year data. The user can choose any date range

startdate : 2012-01-01
enddate: 2012-02-01



dummy - date = startdate  date = enddate

Uwe Ligges



startdate: 2012-06-01
enddate=2012-07-01

the columns of dataset= X ,Y ,date


Thanks and re

-
Thanks in Advance
 Arun
--
View this message in context: 
http://r.789695.n4.nabble.com/Creating-a-dummy-variable-tp4542150p4542150.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] Difference between spec.pgram spec.ar

2012-04-09 Thread Bazman76
OK so I neeed to understan better what it it they are trying to measure.

I understood (incorrectly it seems) that they were simply different methods
to get the same result?

--
View this message in context: 
http://r.789695.n4.nabble.com/Difference-between-spec-pgram-spec-ar-tp4541358p4542985.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] Difference between spec.pgram spec.ar

2012-04-09 Thread Uwe Ligges



On 09.04.2012 17:01, Bazman76 wrote:

OK so I neeed to understan better what it it they are trying to measure.

I understood (incorrectly it seems) that they were simply different methods
to get the same result?


Yes. Also note this is a mailing list and you are lucky I was able to 
remember the context. Please always quote the prior thread and do not 
only send to the list.


Uwe Ligges





--
View this message in context: 
http://r.789695.n4.nabble.com/Difference-between-spec-pgram-spec-ar-tp4541358p4542985.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] How to find best parameter values using deSolve n optim() ?

2012-04-09 Thread Thomas Petzoldt

On 4/5/2012 6:32 PM, mhimanshu wrote:

Hi Thomas,

Thank you so much for your suggestion.
I tried your code and it is working fine. Now when I change the values of Y
in yobs I am getting so many warnings.

say,
yobs- data.frame(
  time = 0:7,
  Y = c(0.00, 3.40, 4.60 ,5.80, 5.80, 6.00, 6.00 ),
  Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)


This data set is wrong because time and Z have 8 elements but Y has only 
7. Let's assume the missing Y is 6.0.



So when i fit the model with the same code that you have written, i got the
following warnings:
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
   such that in the machine, T + H = T on the next step
  (H = step size). Solver will continue anyway.
   In above,  R1 =  0.1484502806322D+01   R2 =  0.2264549048113D-16

and I have got so many such warnings.


This means that the step adjustment algorithm was unable to determine an 
appropriate step size (h), possibly because of wrongly negative 
parameter values.



Can you explain me why this is happening??  and Secondly,  I dont understand
why i am getting parameters values in negatives after fitting??


You can restrict the parameter ranges by using the optional arguments 
upper and lower, see ?modFit for details.



Can you
please help me out in this... :)


The reason for this can be manifold, e.g. wrong model specification, 
unrealistic data, inappropriate accuracy settings, inadequate start 
parameters or just unidentifiable parameters, e.g. if some of them 
depend strongly on each other.


The example below showed clearly that the parameters are strongly 
correlated (0.998 or even 1.000 !!!). This means that not all parameters 
can be identified simultaneously because of high interdependency (see my 
comment in my first post), so you should try to reduce the number of 
fitted parameters. Note however, that Ymax needs to be fitted (or 
adjusted) as well.


Fitting nonlinear models can be an art and requires mathematical 
understanding of the applied techniques (differential equations, 
identifiability analysis, numerical optimization), a good understanding 
of the properties of your (differential equation) model -- and some 
feeling and experience.



Thomas

##-

library(deSolve)
library(FME)

## function derivs
derivs - function(time, y, pars) {
  with (as.list(c(y, pars)), {
dy = a1 * Y * (1 - Y/Ymax) - a2 * (Y + 0.001)
dz = a3 * Y - a4 * Z;
return(list(c(dy, dz)))
  })
}

## parameters
pars - c(a1 = 0.9, a2 = 0.7, a3 = 0.06, a4 = 0.02, Ymax=0.8)
#Ymax - c(0.8)
## initial values
y - c(Y = 0.2, Z = 0.1)
## time steps
times - c(seq(0, 10, 0.1))
## numerical solution of ODE
out - ode(y = y, parms = pars, times = times, func = derivs)


yobs - data.frame(
 time = 0:7,
 Y = c(0.00, 3.40, 4.60 ,5.80, 5.80, 6.00, 6.00, 6.00),
 Z = c(0.1, 0.11, 0.119, 0.128, 0.136, 0.145, 0.153, 0.16)
)

modelCost - function(p) {
  out - ode(y = y, func = derivs, parms = p, times = yobs$time)
  return(modCost(out, yobs, weight=mean))
}

## start values for the parameters
startpars - c(a1 = 5, a2 = 0.002, a3 = 0.002, a4 = 0.02, Ymax=7)
plower - c(a1 = 1, a2 = 0.0001, a3 = 0.0001, a4 = 0.0001, Ymax=0.001)
pupper - c(a1 = 10, a2 = 2, a3 = 1, a4 = 1, Ymax=10)


## fit the model; nprint = 1 shows intermediate results
fit - modFit(f = modelCost, p = startpars,
  upper=pupper, lower=plower, control = list(nprint = 1))
summary(fit)

## graphical result
out2 - ode(y = y, parms = startpars, times = times, func = derivs)
out3 - ode(y = y, parms = fit$par, times = times, func = derivs)
plot(out, out2, out3, obs = yobs)

legend(topleft, legend=c(original, startpars, fitted),
  col = 1:3, lty = 1:3)

summary(fit)

#Parameters:
#  Estimate Std. Error t value Pr(|t|)
#a1   5.610e+00  2.913e+03   0.0020.998
#a2   2.000e+00  2.908e+03   0.0010.999
#a3   1.872e-03  3.162e-03   0.5920.566
#a4   1.188e-04  1.224e-01   0.0010.999
#Ymax 8.908e+00  4.615e+03   0.0020.998
#
#Residual standard error: 0.08001 on 11 degrees of freedom
#
#Parameter correlation:
#   a1   a2   a3  a4 Ymax
#a11.0  1.0 -0.09916 -0.1031  1.0
#a21.0  1.0 -0.09917 -0.1032  1.0
#a3   -0.09916 -0.09917  1.0  0.9981 -0.09917
#a4   -0.10315 -0.10315  0.99807  1. -0.10316
#Ymax  1.0  1.0 -0.09917 -0.1032  1.0

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


Re: [R] Difference between spec.pgram spec.ar

2012-04-09 Thread peter dalgaard

On Apr 9, 2012, at 16:55 , Uwe Ligges wrote:

 
 
 On 08.04.2012 20:39, Bazman76 wrote:
 Hi there,
 
 Can someone explain what the difference between spec.pgram and spec.ar is?
 
 I understand that they attempt to do the same thing one using an AR
 estimation of the underlying series to estimate teh sensity the other using
 the FFT. However when applied to teh same data set they seem to be giving
 quite different results?
 
 http://r.789695.n4.nabble.com/file/n4541358/R_example.jpg
 
 
 Clearly the spec.ar() result seems to be smoothed but the results are also
 generally very different only really sharing the peak as the frequencies go
 towards zero.
 
 Can someone please explain why these two functions produce such different
 results on the same data set?
 
 Because they really measure different things? Why do you expect to get the 
 same output in time as well as in frequency domain?


That wasn't the question There is a raw series and two spectra, and it is 
the difference between the latter that is remarkable. Offhand, this is of 
course a bit like comparing a model fit to a nonparametric smoother of ordinary 
measurements: Sometimes the data just do something that the model says that 
they can't do and you get huge discrepancies. So the first reaction must be 
that an autoregressive model is just wrong for these data. My second reaction 
would be that the original series look a bit like they might have an asymmetric 
distribution, so perhaps a log or a square root transformation is warranted. 

That being said, I'm still quite stumped trying to explain the pattern in the 
raw periodogram. Notice that what you are seeing is not so much an initial peak 
as a region in which the spectrum drops effectively to zero. Have the data by 
any chance been subject to some sort of preprocessing that removes 
low-frequency components?

At any rate, this isn't really about R, is it?

-- 
Peter Dalgaard, Professor,
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] Difference between spec.pgram spec.ar

2012-04-09 Thread Bazman76
oops sorry

n 08.04.2012 20:39, Bazman76 wrote: 
 Hi there, 
 
 Can someone explain what the difference between spec.pgram and spec.ar is? 
 
 I understand that they attempt to do the same thing one using an AR 
 estimation of the underlying series to estimate teh sensity the other
 using 
 the FFT. However when applied to teh same data set they seem to be giving 
 quite different results? 
 
 http://r.789695.n4.nabble.com/file/n4541358/R_example.jpg
 
 
 Clearly the spec.ar() result seems to be smoothed but the results are also 
 generally very different only really sharing the peak as the frequencies
 go 
 towards zero. 
 
 Can someone please explain why these two functions produce such different 
 results on the same data set? 
« [hide part of quote]

*Because they really measure different things? Why do you expect to get 
the same output in time as well as in frequency domain? *

Uwe Ligges 


 
 code shown below: 
 
   library(waveslim) 
 
 vols=read.csv(file=C:/Users/ocuk/My Documents/Abs Vol.csv, header=TRUE, 
 sep=,) 
 x-ts(vols[,1]) 
 #x 
 
 ## LA(8) 
 vol.la8- mra(x, la8, 4, modwt) 
 names(vol.la8)- c(d1, d2, d3, d4, s4) 
 ## plot multiresolution analysis of IBM data 
 #par(mfcol=c(6,1), pty=m, mar=c(5-2,4,4-2,2)) 
 par(mfcol=c(3,1), pty=m, mar=c(5-2,4,4-2,2)) 
 plot.ts(x, axes=F, ylab=, main=(abs rtns)) 
 #for(i in 1:5) 
 # plot.ts(vol.la8[[i]], axes=F, ylab=names(vol.la8)[i]) 
 #axis(side=1, at=seq(0,1600,by=100), 
 # labels=c(0,,200,,400,,600,,800,,1000,,1200,,1400,,1600)) 
 
 spectrum(vol.la8[[1]]) 
 specx- spec.ar(x) 
 
 -- 
 [hidden email] mailing list 
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code. 
« [hide part of quote]

 


*OK so I neeed to understan better what it it they are trying to measure. 

I understood (incorrectly it seems) that they were simply different methods
to get the same result?  *

--
View this message in context: 
http://r.789695.n4.nabble.com/Difference-between-spec-pgram-spec-ar-tp4541358p4543140.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 with Book example of Matrix application

2012-04-09 Thread James Lenihan
I found this example in an Introductory R book in the chapter on Matrices and 
Arrays

The array is 
 m
 [,1] [,2] [,3] [,4] [,5]
[1,]0   12   138   20
[2,]   120   15   28   88
[3,]   13   15069
[4,]8   2860   33
[5,]   20   889   330

The code is  

#returns the minimum value of d[i,j], i !=j and the row attaining 
#the minimum, for square symmetric d; no special policy on ties  
mind  - function(d) {
 n  - nrow(d)
 # add a column to identify row number for apply()
 dd  - cbind(d,1:n)
 wmins - apply(dd[-n,],1,imin)
 # wins will be 2 X n, 1st row being indices and 2nd being values
 i - which.min(wmins[2,])
 j - wmins[1,i]
 return(c(d[i,j],i,j))
}

#this finds the location, value of the minimum in a row x
imin  - function(x) {
   lx  - length(x)
i  - x[lx] # original row number
j -which.min(x[(x+1):(lx-1)]) 
k - i+j
return(c(k,x[k]))
}
The result s/b
[1] 6 3 4

I am getting:

Warning messages:
1: In (x + 1):(lx - 1) :
  numerical expression has 6 elements: only the first used
2: In (x + 1):(lx - 1) :
  numerical expression has 6 elements: only the first used
3: In (x + 1):(lx - 1) :
  numerical expression has 6 elements: only the first used
4: In (x + 1):(lx - 1) :
  numerical expression has 6 elements: only the first used

I have check my typing a number of times. 

Does anyone see an error?

Thanks,

Jim L.
[[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] Help with Book example of Matrix application

2012-04-09 Thread Duncan Murdoch

On 09/04/2012 12:32 PM, James Lenihan wrote:

I found this example in an Introductory R book in the chapter on Matrices and
Arrays


You should probably check with the author of the book (there might be an 
errata page posted somewhere), but it looks like a typo, in your code or 
the original:


The array is
  m
  [,1] [,2] [,3] [,4] [,5]
[1,]0   12   138   20
[2,]   120   15   28   88
[3,]   13   15069
[4,]8   2860   33
[5,]   20   889   330

The code is

#returns the minimum value of d[i,j], i !=j and the row attaining
#the minimum, for square symmetric d; no special policy on ties
mind- function(d) {
  n- nrow(d)
  # add a column to identify row number for apply()
  dd- cbind(d,1:n)
  wmins- apply(dd[-n,],1,imin)
  # wins will be 2 X n, 1st row being indices and 2nd being values
  i- which.min(wmins[2,])
  j- wmins[1,i]
  return(c(d[i,j],i,j))
}

#this finds the location, value of the minimum in a row x
imin- function(x) {
lx- length(x)
 i- x[lx] # original row number
 j-which.min(x[(x+1):(lx-1)])


That line doesn't make sense.  Putting i in place of the innermost x 
would make more sense, but wouldn't work in all cases:  if i == lx, 
you'll get a garbage answer.


Duncan Murdoch

 k- i+j
 return(c(k,x[k]))
}
The result s/b
[1] 6 3 4

I am getting:

Warning messages:
1: In (x + 1):(lx - 1) :
   numerical expression has 6 elements: only the first used
2: In (x + 1):(lx - 1) :
   numerical expression has 6 elements: only the first used
3: In (x + 1):(lx - 1) :
   numerical expression has 6 elements: only the first used
4: In (x + 1):(lx - 1) :
   numerical expression has 6 elements: only the first used

I have check my typing a number of times.

Does anyone see an error?

Thanks,

Jim L.
[[alternative HTML version deleted]]

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


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


Re: [R] Help with Book example of Matrix application

2012-04-09 Thread David Winsemius


On Apr 9, 2012, at 12:32 PM, James Lenihan wrote:

I found this example in an Introductory R book in the chapter on  
Matrices and

Arrays

The array is

m

[,1] [,2] [,3] [,4] [,5]
[1,]0   12   138   20
[2,]   120   15   28   88
[3,]   13   15069
[4,]8   2860   33
[5,]   20   889   330

The code is

#returns the minimum value of d[i,j], i !=j and the row attaining
#the minimum, for square symmetric d; no special policy on ties


Would you consider a more compact operation? That code seems very unR- 
ish.


 min( m[row(m) != col(m)] )
#[1] 6
 which(m ==  min( m[row(m) != col(m)] ), arr.ind=TRUE)
#
 row col
[1,]   4   3
[2,]   3   4
---
 c(min= min( m[row(m) != col(m)] ),
   which(m ==  min( m[row(m) != col(m)] ), arr.ind=TRUE)[1,])

min row col
  6   4   3


I suppose the question could be how to debug that code, in which case  
you should be manually stepping through it to see what the values  
are during the implicit loop.



mind  - function(d) {
n  - nrow(d)
# add a column to identify row number for apply()
dd  - cbind(d,1:n)
wmins - apply(dd[-n,],1,imin)
# wins will be 2 X n, 1st row being indices and 2nd being values
i - which.min(wmins[2,])
j - wmins[1,i]
return(c(d[i,j],i,j))
}

#this finds the location, value of the minimum in a row x
imin  - function(x) {
  lx  - length(x)
   i  - x[lx] # original row number
   j -which.min(x[(x+1):(lx-1)])
   k - i+j
   return(c(k,x[k]))
}
The result s/b
[1] 6 3 4

I am getting:

Warning messages:
1: In (x + 1):(lx - 1) :
 numerical expression has 6 elements: only the first used
2: In (x + 1):(lx - 1) :
 numerical expression has 6 elements: only the first used
3: In (x + 1):(lx - 1) :
 numerical expression has 6 elements: only the first used
4: In (x + 1):(lx - 1) :
 numerical expression has 6 elements: only the first used

I have check my typing a number of times.

Does anyone see an error?

Thanks,

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


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] Appropriate method for sharing data across functions

2012-04-09 Thread Hadley Wickham
 Make OPCON an environment and pass it into the functions that may read it or 
 alter it.  There
 is no real need to pass it out, since environments are changed in-place 
 (unlike lists).  E.g.,
   x - list2env(list(one=1, two=ii, three=3))
   x
  environment: 0x03110890
   objects(x)
  [1] one   three two
   x[[two]]
  [1] ii
   with(x, three+one)
  [1] 4
   f - function(z, env) { env[[newZ]] - z ; sqrt(z) }
   f(10, x)
  [1] 3.162278
   x[[newZ]] # put there by f()
  [1] 10

The advantage of using a reference class is you get an object with
defined behaviour, and something that you can document more easily.

Hadley

-- 
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] Appropriate method for sharing data across functions

2012-04-09 Thread Duncan Murdoch

On 05/04/2012 4:20 PM, John C Nash wrote:

In trying to streamline various optimization functions, I would like to have a 
scratch pad
of working data that is shared across a number of functions. These can be 
called from
different levels within some wrapper functions for maximum likelihood and other 
such
computations. I'm sure there are other applications that could benefit from 
this.

Below are two approaches. One uses the- assignment to a structure I call 
OPCON. The
other attempts to create an environment with this name, but fails. Though I 
have looked at
a number of references, I have so far not found an adequate description of how 
to specify
where the OPCON environment is located. (Both the green and blue books do not 
cover this
topic, at least not under environment in the index.)

Is there a recommended approach to this?



The one I would use is to create all of the functions that need access 
to the scratch pad as local functions within another.  For example,


makethem - function() {
  MAXIMIZE - TRUE
  PARSCALE - rep(1, npar)
  KFN - 0

  ... etc ...

  add1 - function(){
KFN - 1 + KFN
  }

  list(add1 = add1, ... other functions here...)
}

A single call to makethem() will return a list of functions which all 
have access to MAXIMIZE, PARSCALE, etc.  If you like, you can pull them 
out of the list to be called independently, e.g.


fns - makethem()
add1 - fns$add1

and they'll still have access.  (If this is a one-off need, you can use 
local( ) around the definitions, but I prefer the wrapper-function style.



I realize I could use argument lists, but they
get long and tedious with the number of items I may need to pass, though 
passing the OPCON
structure in and out might be the proper way. An onAttach() approach was 
suggested by Paul
Gilbert and tried, but it has so far not succeeded and, unfortunately, does not 
seem to be
usable from source() i.e., cannot be interpreted but must be built first.

JN

Example using-

rm(list=ls())
optstart-function(npar){ # create structure for optimization computations
 # npar is number of parameters ?? test??
 OPCON-list(MAXIMIZE=TRUE, PARSCALE=rep(1,npar), FNSCALE=1,
 KFN=0, KGR=0, KHESS=0)
 # may be other stuff
 ls(OPCON)
}


The code above creates OPCON as a global variable; that's dangerous, if 
you have multiple different optimizers potentially being used.


add1-function(){
 OPCON$KFN-1+OPCON$KFN
 test-OPCON$KFN
}

OPCON-list(MAXIMIZE=TRUE, PARSCALE=rep(1,4), FNSCALE=1,
 KFN=0, KGR=0, KHESS=0)
ls(OPCON)
print(add1())
print(add1())
print(ls.str())

rm(OPCON) # Try to remove the scratchpad
print(ls())

tmp-readline(Now try from within a function)
setup-optstart(4) # Need to sort out how to set this up appropriately
cat(setup =)
print(setup)

print(add1())
print(add1())

rm(OPCON) # Try to remove the scratchpad

==
Example (failing) using new.env:

rm(list=ls())
optstart-function(npar){ # create structure for optimization computations
 # npar is number of parameters ?? test??
 OPCON-new.env(parent=globalenv())
 OPCON-list(MAXIMIZE=TRUE, PARSCALE=rep(1,npar), FNSCALE=1,
 KFN=0, KGR=0, KHESS=0)


This failed because the second assignment wiped out the first:  you 
created a new environment, then threw it away when you assigned a list 
to OPCON.  The other problem with this approach is that you never saved 
OPCON anywhere, so when optstart() was done, it would disappear.


Hadley suggested using a reference class, and Bill suggested using an 
environment explicitly; those are both possible alternatives.  My 
version creates an environment but never names it; it's the environment 
of the add1() function and other functions defined within makethem().  
I'd say the more formal versions they suggested would be better if you 
had a longer (or open-ended)  list of possible functions you wanted to 
work with, or if you wanted to pass your OPCON to the user to manipulate.


Duncan Murdoch


 # may be other stuff
 ls(OPCON)
}

add1-function(){
 OPCON$KFN-1+OPCON$KFN
 test-OPCON$KFN
}

OPCON-new.env(parent=globalenv())
OPCON-list(MAXIMIZE=TRUE, PARSCALE=rep(1,4), FNSCALE=1,
 KFN=0, KGR=0, KHESS=0)
ls(OPCON)
print(add1())
print(add1())
print(ls.str())

rm(OPCON) # Try to remove the scratchpad
print(ls())

tmp-readline(Now try from within a function)
setup-optstart(4) # Need to sort out how to set this up appropriately
cat(setup =)
print(setup)

print(add1())
print(add1())

rm(OPCON) # Try to remove the scratchpad

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

[R] Overall model significance for poisson GLM

2012-04-09 Thread Pat Wilkins
Greetings,

I am running glm models for species counts using a poisson link function.
Normal summary functions for this provide summary statistics in the form of
the deviance, AIC, and p-values for individual predictors.  I would like to
obtain the p-value for the overall model.  So far, I have been using an
analysis of deviance table to check a model against the null model with the
intercept as the only predictor.

Any advice on other methods to obtain the proper p-value would be
appreciated.

Thanks,

Pat

[[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] Difference between spec.pgram spec.ar

2012-04-09 Thread Bazman76
Yes I agree, there may be something pathalogical in the way at least one of
the models handles the data.  That's why I was trying to get a better handle
on how the two functions spec.prgm() and spec.ar() work.

The data has been processed by a wavelet analysis, so what you are seeing as
the raw data is in fact the level1 details from the wavelet anlaysis.

In theory this should only have high frequency components, that was why I am
so surpirsed to see such strong components at low frequencies.

That is not a R quesiton per se, but surely how spec.prgm() and spec.ar()
work is?



--
View this message in context: 
http://r.789695.n4.nabble.com/Difference-between-spec-pgram-spec-ar-tp4541358p4543191.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] binned tabulation

2012-04-09 Thread Delia Shelton



Hi,

I am attempting to tabulate binned data. The '1' represents the appearance of 
the focal mouse pup, and '2' represents the disappearance of the focal mouse 
pup. The code written below is intended to calculate the total time spent 
appeared out of 3600s. For Sample 1, both the hand calculation and R code yield 
the same result, 50. A problem seems to occur when '1' is the last entry. For 
Sample 2, the total time appeared is 53  (hand calculation), however, using the 
R code below yields 55. If you have any suggestions for solving the problem, 
please let me know. 

Thank you in advance for any assistance you may provide.

Delia


Sample 1

0.0 1
40 2
45 1
55 2


Sample 2

0.0 1
40 2
45 1
55 2
57 1



name = read.table(file.choose(),header=F) # opening a data file
colnames(name)-c(Time, Behavior)
name = data.frame(name$Behavior, name$Time)
colnames(name)-c(Behavior, Time)    
name-name[name$Time  3600, ];

## run file 

x-seq(0,3600, by = 60) # total time partition by time which is 60

if (tail(name$Behavior, 1) == 1) {name-rbind(name, c(4, 3600))} else
{name-rbind(name, c(1, 3600))} 

if (nrow(name) %% 2 != 0)
 {name -name[-c(nrow(name)), -c(nrow(name))]}

q-NULL
for (y in (1: nrow(name)))
{
    if (y %% 2 != 0)
    q-c(q, (c(name$Time[y]:name$Time[y +1])))
}

b-table(cut(q,x))

sum(b)

[[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] Overall model significance for poisson GLM

2012-04-09 Thread Bert Gunter
Not really an R question. Post on a statistics site like
http://stats.stackexchange.com/

-- Bert

On Mon, Apr 9, 2012 at 9:51 AM, Pat Wilkins pwilk...@illinois.edu wrote:
 Greetings,

 I am running glm models for species counts using a poisson link function.
 Normal summary functions for this provide summary statistics in the form of
 the deviance, AIC, and p-values for individual predictors.  I would like to
 obtain the p-value for the overall model.  So far, I have been using an
 analysis of deviance table to check a model against the null model with the
 intercept as the only predictor.

 Any advice on other methods to obtain the proper p-value would be
 appreciated.

 Thanks,

 Pat

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



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

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


Re: [R] Difference between spec.pgram spec.ar

2012-04-09 Thread Bert Gunter
On Mon, Apr 9, 2012 at 9:27 AM, Bazman76 h_a_patie...@hotmail.com wrote:
 Yes I agree, there may be something pathalogical in the way at least one of
 the models handles the data.  That's why I was trying to get a better handle
 on how the two functions spec.prgm() and spec.ar() work.

 The data has been processed by a wavelet analysis, so what you are seeing as
 the raw data is in fact the level1 details from the wavelet anlaysis.

 In theory this should only have high frequency components, that was why I am
 so surpirsed to see such strong components at low frequencies.

 That is not a R quesiton per se, but surely how spec.prgm() and spec.ar()
 work is?

Not necessarily, if they are e.g. C or Fortran programs merely called
by R. Indeed, even if written in R, if the algorithms are the issue,
then that is essentially a statistics/numerical analysis matter, not
an R programming one.

-- Bert




 --
 View this message in context: 
 http://r.789695.n4.nabble.com/Difference-between-spec-pgram-spec-ar-tp4541358p4543191.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.



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

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


[R] Listing the contents of an FTP directory via R?

2012-04-09 Thread Jonathan Greenberg
R-helpers:

I'd like to be able to store all the file information from an ftp site
(e.g. file and foldernames) through an R command.  Any ideas how to do
this?  Here's an example site to use:

ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005

--j

-- 
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html

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


Re: [R] Listing the contents of an FTP directory via R?

2012-04-09 Thread steven mosher
A couple of ways.

using Rcurl   you can use the  curlOption of dirlistonly.

otherwise you can read the page and parse.  I've got some code around here
to do that.

Steve

On Mon, Apr 9, 2012 at 11:27 AM, Jonathan Greenberg j...@illinois.eduwrote:

 R-helpers:

 I'd like to be able to store all the file information from an ftp site
 (e.g. file and foldernames) through an R command.  Any ideas how to do
 this?  Here's an example site to use:

 ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005

 --j

 --
 Jonathan A. Greenberg, PhD
 Assistant Professor
 Department of Geography and Geographic Information Science
 University of Illinois at Urbana-Champaign
 607 South Mathews Avenue, MC 150
 Urbana, IL 61801
 Phone: 415-763-5476
 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
 http://www.geog.illinois.edu/people/JonathanGreenberg.html

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


[[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] Difference between spec.pgram spec.ar

2012-04-09 Thread Prof Brian Ripley

On 09/04/2012 18:52, Bert Gunter wrote:

On Mon, Apr 9, 2012 at 9:27 AM, Bazman76h_a_patie...@hotmail.com  wrote:

Yes I agree, there may be something pathalogical in the way at least one of
the models handles the data.  That's why I was trying to get a better handle
on how the two functions spec.prgm() and spec.ar() work.

The data has been processed by a wavelet analysis, so what you are seeing as
the raw data is in fact the level1 details from the wavelet anlaysis.

In theory this should only have high frequency components, that was why I am
so surpirsed to see such strong components at low frequencies.

That is not a R quesiton per se, but surely how spec.prgm() and spec.ar()
work is?


Not necessarily, if they are e.g. C or Fortran programs merely called
by R. Indeed, even if written in R, if the algorithms are the issue,
then that is essentially a statistics/numerical analysis matter, not
an R programming one.


Which is why the help pages often (and do here) have definitive 
references.  So the best course of action is to start following up the 
references.


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

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


Re: [R] Appropriate method for sharing data across functions

2012-04-09 Thread John C Nash

Thanks to Hadley, William and Duncan for suggestions. I'm currently 
implementing a
solution that is close to that of William and Duncan (and learning more about 
environments
in the process). I suspect the reference classes are possibly a more reliable 
long term
solution. I'll plead laziness until I find the other approach won't work.

For anyone interested in these issues, over the next few weeks the outcome of my
investigations and trials should be on the SVN repository of R-forge under 
Optimization
and Solving Tools, https://r-forge.r-project.org/R/?group_id=395, in particular 
within the
packages optfntools and optimx. However, it will be late April before things 
stabilize.

John Nash

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


Re: [R] Listing the contents of an FTP directory via R?

2012-04-09 Thread Jonathan Greenberg
Steven:

Thanks -- I seem to be running into the problem with the link I sent along:

 getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,dirlistonly 
 = TRUE)
Error in function (type, msg, asError = TRUE)  : RETR response: 550

I'm wondering if it might be a passive ftp issue, but:
 getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,ftp.use.epsv=TRUE,
  dirlistonly = TRUE)
Error in function (type, msg, asError = TRUE)  :
  FTP response reading failed

Does not seem to work...  Thoughts?

--j

On Mon, Apr 9, 2012 at 1:32 PM, steven mosher mosherste...@gmail.com wrote:
 A couple of ways.

 using Rcurl   you can use the  curlOption of dirlistonly.

 otherwise you can read the page and parse.  I've got some code around here
 to do that.

 Steve

 On Mon, Apr 9, 2012 at 11:27 AM, Jonathan Greenberg j...@illinois.edu
 wrote:

 R-helpers:

 I'd like to be able to store all the file information from an ftp site
 (e.g. file and foldernames) through an R command.  Any ideas how to do
 this?  Here's an example site to use:

 ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005

 --j

 --
 Jonathan A. Greenberg, PhD
 Assistant Professor
 Department of Geography and Geographic Information Science
 University of Illinois at Urbana-Champaign
 607 South Mathews Avenue, MC 150
 Urbana, IL 61801
 Phone: 415-763-5476
 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
 http://www.geog.illinois.edu/people/JonathanGreenberg.html

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





-- 
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html

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


Re: [R] Listing the contents of an FTP directory via R?

2012-04-09 Thread steven mosher
Ya I hit the same error with my code that reads directories.
I'll try some other stuff . I think I hit this error before I  with usgs.


On Mon, Apr 9, 2012 at 11:40 AM, Jonathan Greenberg j...@illinois.eduwrote:

 Steven:

 Thanks -- I seem to be running into the problem with the link I sent along:

  getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,dirlistonly
 = TRUE)
 Error in function (type, msg, asError = TRUE)  : RETR response: 550

 I'm wondering if it might be a passive ftp issue, but:
  getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,ftp.use.epsv=TRUE,
 dirlistonly = TRUE)
 Error in function (type, msg, asError = TRUE)  :
  FTP response reading failed

 Does not seem to work...  Thoughts?

 --j

 On Mon, Apr 9, 2012 at 1:32 PM, steven mosher mosherste...@gmail.com
 wrote:
  A couple of ways.
 
  using Rcurl   you can use the  curlOption of dirlistonly.
 
  otherwise you can read the page and parse.  I've got some code around
 here
  to do that.
 
  Steve
 
  On Mon, Apr 9, 2012 at 11:27 AM, Jonathan Greenberg j...@illinois.edu
  wrote:
 
  R-helpers:
 
  I'd like to be able to store all the file information from an ftp site
  (e.g. file and foldernames) through an R command.  Any ideas how to do
  this?  Here's an example site to use:
 
  ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005
 
  --j
 
  --
  Jonathan A. Greenberg, PhD
  Assistant Professor
  Department of Geography and Geographic Information Science
  University of Illinois at Urbana-Champaign
  607 South Mathews Avenue, MC 150
  Urbana, IL 61801
  Phone: 415-763-5476
  AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
  http://www.geog.illinois.edu/people/JonathanGreenberg.html
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 



 --
 Jonathan A. Greenberg, PhD
 Assistant Professor
 Department of Geography and Geographic Information Science
 University of Illinois at Urbana-Champaign
 607 South Mathews Avenue, MC 150
 Urbana, IL 61801
 Phone: 415-763-5476
 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
 http://www.geog.illinois.edu/people/JonathanGreenberg.html


[[alternative HTML version deleted]]

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


Re: [R] Problems with CRAN package Ryacas

2012-04-09 Thread Martin Maechler
Apropos:

I don't have the problems, the OP had, but on my ubuntu
notebook, Ryacas does not return expressions (just the strings),
and hence

as.expression( yacas-result )

always gives NULL  and e.g. the   demo(Ryacas-Function)
also fails:

  yacas(expression(deriv(BurrCDF(x,c,k
 k*c*x^(c-1)*(x^c+1)^(-(k+1));
  yy - yacas(expression(deriv(BurrCDF(x,c,k
  yy
 k*c*x^(c-1)*(x^c+1)^(-(k+1));
  str(yy)
 List of 2
  $  : NULL
  $ YacasForm: chr k*c*x^(c-1)*(x^c+1)^(-(k+1));
  - attr(*, class)= chr yacas
  

I have 

- yacas 1.2.2  (standard ubuntu package)
- re-installed Ryacas  (manually from source, just now, 
under R 2.15.0 patched)

 sessionInfo()
R version 2.15.0 Patched (2012-04-09 r58947)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=de_CH.UTF-8   LC_NUMERIC=C   LC_TIME=en_US.UTF-8  
 
 [4] LC_COLLATE=de_CH.UTF-8 LC_MONETARY=en_US.UTF-8
LC_MESSAGES=de_CH.UTF-8   
 [7] LC_PAPER=C LC_NAME=C  LC_ADDRESS=C 
 
[10] LC_TELEPHONE=C LC_MEASUREMENT=de_CH.UTF-8 LC_IDENTIFICATION=C  
 

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

other attached packages:
[1] Ryacas_0.2-11  XML_3.9-4  fortunes_1.5-0 sfsmisc_1.0-20

loaded via a namespace (and not attached):
[1] compiler_2.15.0 tools_2.15.0   
 

--

On our Fedora computers at work (with an *older* version of yacas,
installed manually from sources there),
things work fine :

  str(yacas(expression(Factor(x^2-1
 List of 2
  $ text  :  expression((x + 1) * (x - 1))
  $ OMForm: chr [1:15] OMOBJ   OMA OMS cd=\arith1\ 
name=\times\/ OMA ...
  - attr(*, class)= chr yacas
  as.expression(yacas(expression(Factor(x^2-1
 expression((x + 1) * (x - 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] Listing the contents of an FTP directory via R?

2012-04-09 Thread steven mosher
Its kinda gross, but you can just   download.file() and then parse the
result

!DOCTYPE HTML PUBLIC -//IETF//DTD HTML//EN

HTML

HEAD

TITLEFTP directory /MOTA/MCD15A3.005 at e4ftl01.cr.usgs.gov/TITLE

/HEAD

BODY

H2 ID=WinINetFtpDirectoryFTP directory /MOTA/MCD15A3.005 at
e4ftl01.cr.usgs.gov/H2

HR

H4PRE

##

##



 *** WARNING TO USERS OF THIS SYSTEM ***



 THIS COMPUTER SYSTEM, INCLUDING ALL RELATED EQUIPMENT, NETWORKS, AND
NETWORK

 DEVICES (INCLUDING INTERNET ACCESS), IS PROVIDED BY THE DEPARTMENT OF THE

 INTERIOR (DOI) IN ACCORDANCE WITH THE AGENCY POLICY FOR OFFICIAL USE AND

 LIMITED PERSONAL USE.

 ALL AGENCY COMPUTER SYSTEMS MAY BE MONITORED FOR ALL LAWFUL PURPOSES,
INCLUDING

 BUT NOT LIMITED TO, ENSURING THAT USE IS AUTHORIZED, FOR MANAGEMENT OF THE

 SYSTEM, TO FACILITATE PROTECTION AGAINST UNAUTHORIZED ACCESS, AND TO VERIFY

 SECURITY PROCEDURES, SURVIVABILITY AND OPERATIONAL SECURITY. ANY
INFORMATION ON

 THIS COMPUTER SYSTEM MAY BE EXAMINED, RECORDED, COPIED AND USED FOR
AUTHORIZED

 PURPOSES AT ANY TIME.

 ALL INFORMATION, INCLUDING PERSONAL INFORMATION, PLACED OR SENT OVER THIS
SYSTEM

 MAY BE MONITORED, AND USERS OF THIS SYSTEM ARE REMINDED THAT SUCH
MONITORING

 DOES OCCUR. THEREFORE, THERE SHOULD BE NO EXPECTATION OF PRIVACY WITH
RESPECT

 TO USE OF THIS SYSTEM.

 BY LOGGING INTO THIS AGENCY COMPUTER SYSTEM, YOU ACKNOWLEDGE AND CONSENT TO

 THE MONITORING OF THIS SYSTEM. EVIDENCE OF YOUR USE, AUTHORIZED OR
UNAUTHORIZED,

 COLLECTED DURING MONITORING MAY BE USED FOR CIVIL, CRIMINAL,
ADMINISTRATIVE,

 OR OTHER ADVERSE ACTION. UNAUTHORIZED OR ILLEGAL USE MAY SUBJECT YOU TO

 PROSECUTION.



##

##



/PRE/H4

HR

A HREF=..Up to higher level directory/ABRPRE

10/18/2011 05:23PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.04/B2002.07.04/B/A

10/18/2011 05:08PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.08/B2002.07.08/B/A

10/18/2011 05:05PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.12/B2002.07.12/B/A

10/18/2011 06:00PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.16/B2002.07.16/B/A

10/18/2011 06:00PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.20/B2002.07.20/B/A

10/18/2011 06:00PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.24/B2002.07.24/B/A

10/18/2011 06:01PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.07.28/B2002.07.28/B/A

10/18/2011 06:01PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.05/B2002.08.05/B/A

10/18/2011 06:01PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.09/B2002.08.09/B/A

10/18/2011 05:55PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.13/B2002.08.13/B/A

10/18/2011 05:54PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.17/B2002.08.17/B/A

10/18/2011 05:37PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.21/B2002.08.21/B/A

10/18/2011 04:26PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.25/B2002.08.25/B/A

10/18/2011 05:19PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.08.29/B2002.08.29/B/A

10/18/2011 05:19PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.02/B2002.09.02/B/A

10/18/2011 05:14PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.06/B2002.09.06/B/A

10/18/2011 03:39PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.10/B2002.09.10/B/A

10/18/2011 05:23PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.14/B2002.09.14/B/A

10/18/2011 05:25PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.18/B2002.09.18/B/A

10/18/2011 05:25PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.22/B2002.09.22/B/A

10/18/2011 05:21PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.26/B2002.09.26/B/A

10/18/2011 05:59PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.09.30/B2002.09.30/B/A

10/18/2011 06:00PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.04/B2002.10.04/B/A

10/18/2011 04:37PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.08/B2002.10.08/B/A

10/18/2011 04:03PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.12/B2002.10.12/B/A

10/18/2011 06:00PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.16/B2002.10.16/B/A

10/18/2011 05:36PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.20/B2002.10.20/B/A

10/18/2011 03:54PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.24/B2002.10.24/B/A

10/18/2011 03:54PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.10.28/B2002.10.28/B/A

10/18/2011 03:50PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.11.01/B2002.11.01/B/A

10/18/2011 01:25PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.11.05/B2002.11.05/B/A

10/18/2011 01:27PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.11.09/B2002.11.09/B/A

10/18/2011 01:29PM  Directory A
HREF=/MOTA/MCD15A3.005/2002.11.13/B2002.11.13/B/A

10/18/2011 01:31PM  Directory A

Re: [R] Listing the contents of an FTP directory via R?

2012-04-09 Thread Jonathan Greenberg
It seems to be choking on NLST:

require(RCurl)
getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,ftp.use.epsv=TRUE,
dirlistonly = TRUE)

...

 230 Guest login ok, access restrictions apply.
 PWD
 257 / is current directory.
* Entry path is '/'
 CWD MOTA
 250 CWD command successful.
 CWD MCD15A3.005
 250 CWD command successful.
 EPSV
* Connect data stream passively
 500 'EPSV': command not understood.
* disabling EPSV usage
 PASV
 227 Entering Passive Mode (152,61,133,5,177,17)
*   Trying 152.61.133.5... * connected
* Connecting to 152.61.133.5 (152.61.133.5) port 45329
 TYPE A
 200 Type set to A.
 NLST
 550 No files found.
* RETR response: 550
* Remembering we are in dir MOTA/MCD15A3.005/
* Connection #0 to host e4ftl01.cr.usgs.gov left intact
Error in function (type, msg, asError = TRUE)  : RETR response: 550

Thanks for looking into this!

--j

On Mon, Apr 9, 2012 at 1:45 PM, steven mosher mosherste...@gmail.com wrote:
 Ya I hit the same error with my code that reads directories.
 I'll try some other stuff . I think I hit this error before I  with usgs.


 On Mon, Apr 9, 2012 at 11:40 AM, Jonathan Greenberg j...@illinois.edu
 wrote:

 Steven:

 Thanks -- I seem to be running into the problem with the link I sent
 along:

 
  getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,dirlistonly
  = TRUE)
 Error in function (type, msg, asError = TRUE)  : RETR response: 550

 I'm wondering if it might be a passive ftp issue, but:
 
  getURL(ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005/,verbose=TRUE,ftp.use.epsv=TRUE,
  dirlistonly = TRUE)
 Error in function (type, msg, asError = TRUE)  :
  FTP response reading failed

 Does not seem to work...  Thoughts?

 --j

 On Mon, Apr 9, 2012 at 1:32 PM, steven mosher mosherste...@gmail.com
 wrote:
  A couple of ways.
 
  using Rcurl   you can use the  curlOption of dirlistonly.
 
  otherwise you can read the page and parse.  I've got some code around
  here
  to do that.
 
  Steve
 
  On Mon, Apr 9, 2012 at 11:27 AM, Jonathan Greenberg j...@illinois.edu
  wrote:
 
  R-helpers:
 
  I'd like to be able to store all the file information from an ftp site
  (e.g. file and foldernames) through an R command.  Any ideas how to do
  this?  Here's an example site to use:
 
  ftp://e4ftl01.cr.usgs.gov/MOTA/MCD15A3.005
 
  --j
 
  --
  Jonathan A. Greenberg, PhD
  Assistant Professor
  Department of Geography and Geographic Information Science
  University of Illinois at Urbana-Champaign
  607 South Mathews Avenue, MC 150
  Urbana, IL 61801
  Phone: 415-763-5476
  AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
  http://www.geog.illinois.edu/people/JonathanGreenberg.html
 
  __
  R-help@r-project.org mailing list
  https://stat.ethz.ch/mailman/listinfo/r-help
  PLEASE do read the posting guide
  http://www.R-project.org/posting-guide.html
  and provide commented, minimal, self-contained, reproducible code.
 
 



 --
 Jonathan A. Greenberg, PhD
 Assistant Professor
 Department of Geography and Geographic Information Science
 University of Illinois at Urbana-Champaign
 607 South Mathews Avenue, MC 150
 Urbana, IL 61801
 Phone: 415-763-5476
 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
 http://www.geog.illinois.edu/people/JonathanGreenberg.html





-- 
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html

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


Re: [R] Problems with CRAN package Ryacas

2012-04-09 Thread Gabor Grothendieck
On Mon, Apr 9, 2012 at 2:47 PM, Martin Maechler
maech...@stat.math.ethz.ch wrote:
 Apropos:

 I don't have the problems, the OP had, but on my ubuntu
 notebook, Ryacas does not return expressions (just the strings),
 and hence

    as.expression( yacas-result )

 always gives NULL  and e.g. the   demo(Ryacas-Function)
 also fails:

   yacas(expression(deriv(BurrCDF(x,c,k
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   yy - yacas(expression(deriv(BurrCDF(x,c,k
   yy
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   str(yy)
  List of 2
  $          : NULL
  $ YacasForm: chr k*c*x^(c-1)*(x^c+1)^(-(k+1));
  - attr(*, class)= chr yacas
  

 I have

 - yacas 1.2.2  (standard ubuntu package)

The latest version of Ryacas on CRAN is 0.2-11
http://cran.r-project.org/package=Ryacas

-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at 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] Problems with CRAN package Ryacas

2012-04-09 Thread Gabor Grothendieck
On Mon, Apr 9, 2012 at 4:02 PM, Gabor Grothendieck
ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 2:47 PM, Martin Maechler
 maech...@stat.math.ethz.ch wrote:
 Apropos:

 I don't have the problems, the OP had, but on my ubuntu
 notebook, Ryacas does not return expressions (just the strings),
 and hence

    as.expression( yacas-result )

 always gives NULL  and e.g. the   demo(Ryacas-Function)
 also fails:

   yacas(expression(deriv(BurrCDF(x,c,k
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   yy - yacas(expression(deriv(BurrCDF(x,c,k
   yy
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   str(yy)
  List of 2
  $          : NULL
  $ YacasForm: chr k*c*x^(c-1)*(x^c+1)^(-(k+1));
  - attr(*, class)= chr yacas
  

 I have

 - yacas 1.2.2  (standard ubuntu package)

 The latest version of Ryacas on CRAN is 0.2-11
 http://cran.r-project.org/package=Ryacas

and the version of yacas supported by Ryacas is:

 yacas(Version())
expression(1.0.63)

The troubleshooting section on the home page that I referred to does
mention this.

-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at 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.


[R] sdev, variance in prcomp

2012-04-09 Thread carol white
Hello,
It might be a trivial question but I just wanted to find out the relationship 
between sdev and proportion of variance generated by prcomp. I got the 
following result from my data set

 PC1  PC2  PC3
Standard deviation 104.89454 15.40910 9.012047
Proportion of Variance   0.52344  0.01130 0.003860
Cumulative Proportion    0.52344  0.53474 0.538600



first, I had thought that the variance must be standard deviation power to 2 
but it seems that variance is sdev/200. Did I misunderstood some thing?

Best,

Carol


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

2012-04-09 Thread bekir . cetintav
Hi!
I study fuzzy regression. Linear Programming(LP) methods are commonly used
for fuzzy linear regression (FLR) because they are simple and easy to
apply. I should make a simulation study with R  and have come a distance.
I use linp
(for instence a part from my codes)
(L-linp(E=NULL,F=NULL,Cost=obj,G=sınır,H=sol,ispos=FALSE))
there is ispose argument which allows all the variables(x) to get
negative values when its false. my problem is there. I want to allow some
of the variables(x)to get negative values and the others must be
nonnegative.(for instence x1,x2,x3 could be negative or positive and
x4,x5,x6 must be nonnegative.) Is there anything to solve that problem?
thanks for any help.
Bekir.

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


Re: [R] sdev, variance in prcomp

2012-04-09 Thread Kevin Wright
It is reporting _proportion_ of the total variance.  Not the variance.

The total variance is the sum of all the squared standard deviations.  Only
the first 3 standard deviations are shown.

Kevin


On Mon, Apr 9, 2012 at 3:31 PM, carol white wht_...@yahoo.com wrote:

 Hello,
 It might be a trivial question but I just wanted to find out the
 relationship between sdev and proportion of variance generated by prcomp. I
 got the following result from my data set

  PC1  PC2  PC3
 Standard deviation 104.89454 15.40910 9.012047
 Proportion of Variance   0.52344  0.01130 0.003860
 Cumulative Proportion0.52344  0.53474 0.538600



 first, I had thought that the variance must be standard deviation power to
 2 but it seems that variance is sdev/200. Did I misunderstood some thing?

 Best,

 Carol


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




-- 
Kevin Wright

[[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 comparison matrix elements

2012-04-09 Thread Kerapi
Hi!, 
I'm really hoping someone out there will be able to help me.

I recently started my MSc dissertation on Population Projection Matrices, which 
has been going well until now. I am trying to set-up a general script that does 
a pairwise comparison of all elements in my matrices. 

So for example, given that I have the following matrix S:

 S

[,1] [,2] [,3]

[1,] 0.0 0.007361183 0.0001634936

[2,] 13.88458 0.353988378 0.00

[3,] 0.0 16.086367098 0.3572819557

I'd like to create a matrix SA for each element in my matrix compared to 
another element like so:

 SA 

S-[1,1]

 [,1][,2][,3]

[1,] [1,1]-[1,1] [1,2]-[1,1] [1,3]-[1,1]

[2,] [2,1]-[1,1] [2,2]-[1,1] [2,3]-[1,1]

[3,] [3,1]-[1,1] [3,2]-[1,1] [3,3]-[1,1]



S-[1,2]

[...]



S-[1,3]

[...]

... and so on

The aim is to be able to prove existing rules and trends across matrix 
dimensions. 

For example, I've been trying to test whether the first line of values 
decreases or remains the same from column 1 to the penultimate column: S[1,j] 
= S[1,lj]

The last thing I tried was:

M-Matlab2R([0 24.724 1377.48;0.029328 0.26 0;0 0.014 0.78]) 

w - abs(Re(eigen(M)$vectors[,1]))

v - abs(Re(eigen(t(M))$vectors[,1]))

w - w/sum(w)

v - v/(t(v)%*%w) 

S - (v%*%t(w))/as.vector(t(v)%*%w)

S



n- length(S)

M.comp - array(S,dim=c(n,n,n,n))

for (i in 1:n) {

for (j in 1:n) {

for (k in 1:n) {

for (l in 1:n) {

poscompa[1,3,1,j-1] - S[i,j]=S[k,l] # if=1, then TRUE and if=0, then FALSE 



M.comp

sum(poscompa[1,3,1,1:n]) # should give me the value for n. I know this doesn't 
test it accurately, but first I need to get the loops to work right.

I'm getting an array of errors such as can't find function poscompa (or 
compa), subscripts out of range or type of closure not valid regardless of 
what I try. Using COMPAR (poscompa and compa) was the last recommendation I 
got, but I'm starting to wonder if there might be other ways to go about this. 
All out-of-the-box ideas I've come up with and tried haven't gotten me much 
farther. I've now practically exhausted my creative thinking, and I'm becoming 
very frustrated. I'd really like to get this script going since my current one 
would make my life hell for large populations (60x60 population matrices). 

Any ideas on how I could move forward? 



Many, many thanks!

Marta

[[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] binned tabulation

2012-04-09 Thread Jean V Adams
Delia,

name - data.frame(Behavior=c(1, 2, 1, 2, 1), Time=c(0, 40, 45, 55, 57))

appear - name$Time[name$Behavior==1]
disappear - name$Time[name$Behavior==2]
if(length(appear)  length(disappear)) disappear - c(disappear, 60)
sum(disappear - appear)

Delia Shelton wrote on 04/09/2012 12:30:23 PM:

 Hi,
 
 I am attempting to tabulate binned data. The '1' represents the 
 appearance of the focal mouse pup, and '2' represents the 
 disappearance of the focal mouse pup. The code written below is 
 intended to calculate the total time spent appeared out of 3600s. 
 For Sample 1, both the hand calculation and R code yield the same 
 result, 50. A problem seems to occur when '1' is the last entry. For
 Sample 2, the total time appeared is 53  (hand calculation), 
 however, using the R code below yields 55. If you have any 
 suggestions for solving the problem, please let me know. 
 
 Thank you in advance for any assistance you may provide.
 
 Delia
 
 
 Sample 1
 
 0.0 1
 40 2
 45 1
 55 2
 
 
 Sample 2
 
 0.0 1
 40 2
 45 1
 55 2
 57 1
 
 
 
 name = read.table(file.choose(),header=F) # opening a data file
 colnames(name)-c(Time, Behavior)
 name = data.frame(name$Behavior, name$Time)
 colnames(name)-c(Behavior, Time)
 name-name[name$Time  3600, ];
 
 ## run file 
 
 x-seq(0,3600, by = 60) # total time partition by time which is 60
 
 if (tail(name$Behavior, 1) == 1) {name-rbind(name, c(4, 3600))} else
 {name-rbind(name, c(1, 3600))} 
 
 if (nrow(name) %% 2 != 0)
  {name -name[-c(nrow(name)), -c(nrow(name))]}
 
 q-NULL
 for (y in (1: nrow(name)))
 {
 if (y %% 2 != 0)
 q-c(q, (c(name$Time[y]:name$Time[y +1])))
 }
 
 b-table(cut(q,x))
 
 sum(b)

[[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] Pairwise comparison matrix elements

2012-04-09 Thread David Winsemius


On Apr 9, 2012, at 3:17 PM, Kerapi wrote:


Hi!,
I'm really hoping someone out there will be able to help me.

I recently started my MSc dissertation on Population Projection  
Matrices, which has been going well until now. I am trying to set-up  
a general script that does a pairwise comparison of all elements in  
my matrices.


So for example, given that I have the following matrix S:


S


   [,1] [,2] [,3]

[1,] 0.0 0.007361183 0.0001634936

[2,] 13.88458 0.353988378 0.00

[3,] 0.0 16.086367098 0.3572819557

I'd like to create a matrix SA for each element in my matrix  
compared to another element like so:



SA


S-[1,1]

[,1][,2][,3]

[1,] [1,1]-[1,1] [1,2]-[1,1] [1,3]-[1,1]

[2,] [2,1]-[1,1] [2,2]-[1,1] [2,3]-[1,1]

[3,] [3,1]-[1,1] [3,2]-[1,1] [3,3]-[1,1]



Try this:

kS - -1*kronecker(c(S),S, -)
dim(kS) - c(3, 3, 3, 3)

The row number is the thrid dimension and the col nimber is the second  
dimension. If you want them to be rearrange to be in the ( r,c) order  
you can use aperm.


 kSr - aperm(kS, c(1, 4, 2, 3))

#The zero entires were used because I could easily tell what a correct  
answer should be.

 kSr[,,3,1]
 [,1] [,2] [,3]
[1,]  0.0  0.007361183 0.0001634936
[2,] 13.88458  0.353988378 0.00
[3,]  0.0 16.086367098 0.3572819557
 kSr[,,1,1]
 [,1] [,2] [,3]
[1,]  0.0  0.007361183 0.0001634936
[2,] 13.88458  0.353988378 0.00
[3,]  0.0 16.086367098 0.3572819557
 kSr[,,2,3]
 [,1] [,2] [,3]
[1,]  0.0  0.007361183 0.0001634936
[2,] 13.88458  0.353988378 0.00
[3,]  0.0 16.086367098 0.3572819557

Test a non-zero entry
 kSr[,,2,2]
   [,1]   [,2] [,3]
[1,] -0.3539884 -0.3466272 -0.353824884
[2,] 13.5305916  0.000 -0.353988378
[3,] -0.3539884 15.7323787  0.003293578
 S - S[2,2]
   [,1]   [,2] [,3]
[1,] -0.3539884 -0.3466272 -0.353824884
[2,] 13.5305916  0.000 -0.353988378
[3,] -0.3539884 15.7323787  0.003293578

 sapply(1:3, function(rw) sapply(1:3,
  function(cl) identical(S-S[rw,cl],  kSr[,,rw,cl] )) )
 [,1] [,2] [,3]
[1,] TRUE TRUE TRUE
[2,] TRUE TRUE TRUE
[3,] TRUE TRUE TRUE

--
David.



S-[1,2]

[...]



S-[1,3]

[...]

... and so on

The aim is to be able to prove existing rules and trends across  
matrix dimensions.


For example, I've been trying to test whether the first line of  
values decreases or remains the same from column 1 to the  
penultimate column: S[1,j] = S[1,lj]


The last thing I tried was:

M-Matlab2R([0 24.724 1377.48;0.029328 0.26 0;0 0.014 0.78])

w - abs(Re(eigen(M)$vectors[,1]))

v - abs(Re(eigen(t(M))$vectors[,1]))

w - w/sum(w)

v - v/(t(v)%*%w)

S - (v%*%t(w))/as.vector(t(v)%*%w)

S



n- length(S)

M.comp - array(S,dim=c(n,n,n,n))

for (i in 1:n) {

for (j in 1:n) {

for (k in 1:n) {

for (l in 1:n) {

poscompa[1,3,1,j-1] - S[i,j]=S[k,l] # if=1, then TRUE and if=0,  
then FALSE




M.comp

sum(poscompa[1,3,1,1:n]) # should give me the value for n. I know  
this doesn't test it accurately, but first I need to get the loops  
to work right.


I'm getting an array of errors such as can't find function poscompa  
(or compa), subscripts out of range or type of closure not  
valid regardless of what I try. Using COMPAR (poscompa and compa)  
was the last recommendation I got, but I'm starting to wonder if  
there might be other ways to go about this. All out-of-the-box ideas  
I've come up with and tried haven't gotten me much farther. I've now  
practically exhausted my creative thinking, and I'm becoming very  
frustrated. I'd really like to get this script going since my  
current one would make my life hell for large populations (60x60  
population matrices).


Any ideas on how I could move forward?



Many, many thanks!

Marta

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


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] For loops

2012-04-09 Thread Christopher Desjardins
Hi,
I am having trouble with syntax for a for loop. Here is what I am trying to
do.

class=c(rep(1,3),rep(2,3),rep(3,3))
out1=rnorm(length(class))
out2=rnorm(length(class))
out3=rnorm(length(class))
data=data.frame(class,out1,out2,out3)

dat.split=split(data,data$class)
  for(i in 1:3){
  sub[i]=dat.split[i]
  }

However, the for loop doesn't work. I want to assign each split to a
different data object. Better yet, how I could assign each class to a
separate object and skip the splitting?

Thanks,
Chris

[[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] For loops

2012-04-09 Thread Jeff Newmiller
class1 - data[1==data$class,]

gets you a subset of the data into a dedicated object, but if you want to 
handle arbitrarily large amounts of data or class values then the list output 
of split is really much better to stay with.

Also, data is a predefined function, so it is not a good idea to define 
objects with that name due to confusing code review.
---
Jeff NewmillerThe .   .  Go Live...
DCN:jdnew...@dcn.davis.ca.usBasics: ##.#.   ##.#.  Live Go...
  Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/BatteriesO.O#.   #.O#.  with
/Software/Embedded Controllers)   .OO#.   .OO#.  rocks...1k
--- 
Sent from my phone. Please excuse my brevity.

Christopher Desjardins cddesjard...@gmail.com wrote:

Hi,
I am having trouble with syntax for a for loop. Here is what I am
trying to
do.

class=c(rep(1,3),rep(2,3),rep(3,3))
out1=rnorm(length(class))
out2=rnorm(length(class))
out3=rnorm(length(class))
data=data.frame(class,out1,out2,out3)

dat.split=split(data,data$class)
  for(i in 1:3){
  sub[i]=dat.split[i]
  }

However, the for loop doesn't work. I want to assign each split to a
different data object. Better yet, how I could assign each class to a
separate object and skip the splitting?

Thanks,
Chris

   [[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] For loops

2012-04-09 Thread David Winsemius


On Apr 9, 2012, at 5:33 PM, Christopher Desjardins wrote:


Hi,
I am having trouble with syntax for a for loop. Here is what I am  
trying to

do.

class=c(rep(1,3),rep(2,3),rep(3,3))
out1=rnorm(length(class))
out2=rnorm(length(class))
out3=rnorm(length(class))
data=data.frame(class,out1,out2,out3)

dat.split=split(data,data$class)
 for(i in 1:3){
 sub[i]=dat.split[i]
 }

However, the for loop doesn't work. I want to assign each split to a
different data object.


Why? What's wrong with leaving them in a list that split() provides.

In that form you can easily use lapply() and start your 12 step  
program away from for-loop addiction.



Better yet, how I could assign each class to a
separate object and skip the splitting?


--

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] Overall model significance for poisson GLM

2012-04-09 Thread Ben Bolker
Pat Wilkins pwilkin2 at illinois.edu writes:

 
 Greetings,
 
 I am running glm models for species counts using a poisson link function.
 Normal summary functions for this provide summary statistics in the form of
 the deviance, AIC, and p-values for individual predictors.  I would like to
 obtain the p-value for the overall model.  So far, I have been using an
 analysis of deviance table to check a model against the null model with the
 intercept as the only predictor.
 
 Any advice on other methods to obtain the proper p-value would be
 appreciated.
 

  What you're doing seems reasonable, although you can also dig
the necessary values out of the summary and compute the p-value
yourself:

counts - c(18,17,15,20,10,20,25,13,12)
outcome - gl(3,1,9)
treatment - gl(3,3)
d.AD - data.frame(treatment, outcome, counts)
glm.D93 - glm(counts ~ outcome + treatment, family=poisson(), data=d.AD)

## as you have been doing
anova(update(glm.D93,.~1),glm.D93,test=Chisq)
## or
sg1 - summary(glm.D93)
devdiff - with(sg1,null.deviance-deviance)
dfdiff - with(sg1,df.null-df.residual)
pchisq(abs(devdiff),df=dfdiff,lower.tail=FALSE)

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

2012-04-09 Thread Hasan Diwan
I have time-series data looking like this:
 dataIn[sample(c(1:nrow(dataIn)), 25),]
 accelerometer_y id data_block_epoch_time
782   0.8424 201300 133179733
1868  0.3432 202386 1331797384000
1828  0.3510 202346 1331797382000
1026  0.2184 201544 1331797342000
1569  0.3432 202087 1331797369000
1453  0.3588 201971 1331797363000
1204  0.3666 201722 1331797351000
1821  0.3588 202339 1331797382000
860   0.8658 201378 1331797333000
910   0.8580 201428 1331797336000
1488  0.3432 202006 1331797365000
578   0.9126 201096 1331797319000
1478  0.3666 201996 1331797364000
1183  0.3588 201701 133179735
29   -0.1716 200547 1331797292000
1540  0.3588 202058 1331797367000
392  -0.1560 200910 133179731
1533  0.3744 202051 1331797367000
1016  0.6318 201534 1331797341000
314  -0.1560 200832 1331797306000
410  -0.1638 200928 1331797311000
769   0.8580 201287 1331797329000
1101  0.3588 201619 1331797346000
403  -0.1638 200921 1331797311000
1794  0.3666 202312 133179738

The id field represents the subsecond value in the timestamp (n
ids/second). What I'd like to do is combine the fields, so that I'm left
with an equally-spaced timeseries with readings. What's the most efficient
way to do this using R? There are an arbitrary number of timestamps per
second. Many thanks! -- H
-- 
Sent from my mobile device
Envoyait de mon portable

[[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] how to add 3d-points to bplot {rms} figure?

2012-04-09 Thread phi771
Hello!
I have created a bplot-figure using this code:

*file - 2dcali_red.ttt
ux-as.matrix(read.table(file, dec = ,))

mode(ux)-'numeric'

vel-ux[,1]
ang-ux[,2] 
x-ux[,3]
y-ux[,4]

dat- data.frame(ang=ang, x=x,y=y) 

require(rms)

ddist2 - datadist(dat)
options(datadist=ddist2)

fitn - lrm(ang ~ rcs(x,4) + rcs(y,4), data=dat)
predi - Predict(fitn, x, y)

bplot(predi,lfun=wireframe, screen = list(z = -40, x = -80), drape=TRUE)*

The file 2dcali_red.ttt  consists of 4 columns can be found here :
http://www.color-space.de/upload/dl/2dcali_red.ttt link 

The code gives me the bplot-figure, which looks like this:
http://r.789695.n4.nabble.com/file/n4544050/Screen_shot_2012-04-10_at_12.04.24_AM.png
 

Now I want to add some 3d data (x,y,z) to this plot (like scatterplot3d),
but i don't know how. I also don't know how to extract the x,y,z mesh as
coordinates from the bplot. Some help would be very appreciated. Thank you!

--
View this message in context: 
http://r.789695.n4.nabble.com/how-to-add-3d-points-to-bplot-rms-figure-tp4544050p4544050.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] Problems with CRAN package Ryacas

2012-04-09 Thread Kjetil Halvorsen
see below!

On Mon, Apr 9, 2012 at 3:04 PM, Gabor Grothendieck
ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 4:02 PM, Gabor Grothendieck
 ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 2:47 PM, Martin Maechler
 maech...@stat.math.ethz.ch wrote:
 Apropos:

 I don't have the problems, the OP had, but on my ubuntu
 notebook, Ryacas does not return expressions (just the strings),
 and hence

    as.expression( yacas-result )

 always gives NULL  and e.g. the   demo(Ryacas-Function)
 also fails:

   yacas(expression(deriv(BurrCDF(x,c,k
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   yy - yacas(expression(deriv(BurrCDF(x,c,k
   yy
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   str(yy)
  List of 2
  $          : NULL
  $ YacasForm: chr k*c*x^(c-1)*(x^c+1)^(-(k+1));
  - attr(*, class)= chr yacas
  

 I have

 - yacas 1.2.2  (standard ubuntu package)

 The latest version of Ryacas on CRAN is 0.2-11
 http://cran.r-project.org/package=Ryacas

 and the version of yacas supported by Ryacas is:

 yacas(Version())
 expression(1.0.63)

 The troubleshooting section on the home page that I referred to does
 mention this.

Thanks!

I deinstalled the synaptic-installed yacas, downloaded yacas_1.0.63.tgz
unpacked it , run
./configure --enable-server
then did make

But make does not terminate, it ends in error! I give the last few
lines of output:

rm -fr .libs/filescanner.la .libs/filescanner.* .libs/filescanner.*
gcc -shared  filescanner.lo plugin.lo   -Wl,-soname -Wl,filescanner.so
-o .libs/filescanner.so
ar cru .libs/filescanner.a  filescanner.o plugin.o
ranlib .libs/filescanner.a
creating filescanner.la
(cd .libs  rm -f filescanner.la  ln -s ../filescanner.la filescanner.la)
make[3]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins/filescanner'
make[3]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
make[3]: Nothing to be done for `all-am'.
make[3]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
Making all in proteus
make[2]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/proteus'
make[2]: Nothing to be done for `all'.
make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/proteus'
Making all in manmake
make[2]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/manmake'
source='manripper.cpp' object='manripper.o' libtool=no \
depfile='.deps/manripper.Po' tmpdepfile='.deps/manripper.TPo' \
depmode=gcc3 /bin/bash ../depcomp \
g++ -DHAVE_CONFIG_H -I. -I. -I.. -g -O2 -Wall  -c -o manripper.o
`test -f 'manripper.cpp' || echo './'`manripper.cpp
manripper.cpp: In function ‘void GetBf(char*, int, FILE*)’:
manripper.cpp:20:21: error: ‘strlen’ was not declared in this scope
manripper.cpp: In function ‘void ProcessFile(char*)’:
manripper.cpp:40:23: error: ‘strlen’ was not declared in this scope
manripper.cpp:43:30: error: ‘strncmp’ was not declared in this scope
manripper.cpp:57:25: error: ‘memset’ was not declared in this scope
manripper.cpp:58:39: error: ‘memcpy’ was not declared in this scope
manripper.cpp:78:27: error: ‘strlen’ was not declared in this scope
manripper.cpp:81:25: error: ‘memset’ was not declared in this scope
manripper.cpp:82:41: error: ‘memcpy’ was not declared in this scope
make[2]: *** [manripper.o] Error 1
make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/manmake'
make[1]: *** [all-recursive] Error 1
make[1]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63'
make: *** [all] Error 2
kjetil@kjetil:~/yacas/yacas-1.0.63$

Kjetil




 --
 Statistics  Software Consulting
 GKX Group, GKX Associates Inc.
 tel: 1-877-GKX-GROUP
 email: ggrothendieck at 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.

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


Re: [R] Gradients in bar charts XXXX

2012-04-09 Thread ilai
On Mon, Apr 9, 2012 at 12:40 PM, Jason Rodriguez
jason.rodrig...@dca.ga.gov wrote:
 Hello, I have a graphics-related question:

 I was wondering if anyone knows of a way to create a bar chart that is 
 colored with a three-part gradient that changes at fixed y-values. Each bar 
 needs to fade green-to-yellow at Y=.10 and from yellow-to-red at Y=.20. Is 
 there an option in a package somewhere that offers an easy way to do this?

?rainbow ?hsv
In R an easy way is an ill-defined term. In the absence of actual data:

bpd - matrix(c(1,seq(0,1,l=64),2,1,seq(0,1,l=64),5,1,seq(0,1,l=64),7),nc=3)
mycols - c('green',rainbow(64,start=0,end=.4)[64:1],'red')
barplot(bpd,col=mycols,border=NA)

Easy enough ?

Cheers



 Attached is a chart I macgyvered together in Excel using a combination of a 
 simple bar chart, fit line, and some drawing tools. I want to avoid doing it 
 this way in the future by finding a way to replicate it in R.

 Any ideas?

 Thanks,

 Jason Michael Rodriguez
 Data Analyst
 State Housing Trust Fund for the Homeless
 Georgia Department of Community Affairs
 Email:  jason.rodrig...@dca.ga.gov


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


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


[R] how to convert seconds to 12 hour time format

2012-04-09 Thread cassie jones
Hello everyone,

I am wondering if there is any routine in R which can convert time given in
'seconds' unit to the 12 hour time format. For example, suppose the data
set looks like

x=c(36885,84000,20)  #x in seconds

I want to get the output as

[1]  11:14:45 AM
[2]  11:20:00 PM
[3] 12:20:00 AM

Does anyone have any idea? Thanks in advance.


Cassie

[[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] Problems with CRAN package Ryacas

2012-04-09 Thread Gabor Grothendieck
We have had Linux users that successfully used Ryacas but if you can't
get yacas to build on your particular system then you could try
running the Windows version under wine.  A Windows binary of yacas is
available so you won't have to build yacas.  I don't know if anyone
has tried that yet but its worth a try.

Another possibility is to try rsympy or rmathpiper (see introductory
paragraphs on Ryacas home page for links).

On Mon, Apr 9, 2012 at 8:00 PM, Kjetil Halvorsen
kjetilbrinchmannhalvor...@gmail.com wrote:
 see below!

 On Mon, Apr 9, 2012 at 3:04 PM, Gabor Grothendieck
 ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 4:02 PM, Gabor Grothendieck
 ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 2:47 PM, Martin Maechler
 maech...@stat.math.ethz.ch wrote:
 Apropos:

 I don't have the problems, the OP had, but on my ubuntu
 notebook, Ryacas does not return expressions (just the strings),
 and hence

    as.expression( yacas-result )

 always gives NULL  and e.g. the   demo(Ryacas-Function)
 also fails:

   yacas(expression(deriv(BurrCDF(x,c,k
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   yy - yacas(expression(deriv(BurrCDF(x,c,k
   yy
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   str(yy)
  List of 2
  $          : NULL
  $ YacasForm: chr k*c*x^(c-1)*(x^c+1)^(-(k+1));
  - attr(*, class)= chr yacas
  

 I have

 - yacas 1.2.2  (standard ubuntu package)

 The latest version of Ryacas on CRAN is 0.2-11
 http://cran.r-project.org/package=Ryacas

 and the version of yacas supported by Ryacas is:

 yacas(Version())
 expression(1.0.63)

 The troubleshooting section on the home page that I referred to does
 mention this.

 Thanks!

 I deinstalled the synaptic-installed yacas, downloaded yacas_1.0.63.tgz
 unpacked it , run
 ./configure --enable-server
 then did make

 But make does not terminate, it ends in error! I give the last few
 lines of output:

 rm -fr .libs/filescanner.la .libs/filescanner.* .libs/filescanner.*
 gcc -shared  filescanner.lo plugin.lo   -Wl,-soname -Wl,filescanner.so
 -o .libs/filescanner.so
 ar cru .libs/filescanner.a  filescanner.o plugin.o
 ranlib .libs/filescanner.a
 creating filescanner.la
 (cd .libs  rm -f filescanner.la  ln -s ../filescanner.la filescanner.la)
 make[3]: Leaving directory 
 `/home/kjetil/yacas/yacas-1.0.63/plugins/filescanner'
 make[3]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
 make[3]: Nothing to be done for `all-am'.
 make[3]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
 make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
 Making all in proteus
 make[2]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/proteus'
 make[2]: Nothing to be done for `all'.
 make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/proteus'
 Making all in manmake
 make[2]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/manmake'
 source='manripper.cpp' object='manripper.o' libtool=no \
        depfile='.deps/manripper.Po' tmpdepfile='.deps/manripper.TPo' \
        depmode=gcc3 /bin/bash ../depcomp \
        g++ -DHAVE_CONFIG_H -I. -I. -I..     -g -O2 -Wall  -c -o manripper.o
 `test -f 'manripper.cpp' || echo './'`manripper.cpp
 manripper.cpp: In function ‘void GetBf(char*, int, FILE*)’:
 manripper.cpp:20:21: error: ‘strlen’ was not declared in this scope
 manripper.cpp: In function ‘void ProcessFile(char*)’:
 manripper.cpp:40:23: error: ‘strlen’ was not declared in this scope
 manripper.cpp:43:30: error: ‘strncmp’ was not declared in this scope
 manripper.cpp:57:25: error: ‘memset’ was not declared in this scope
 manripper.cpp:58:39: error: ‘memcpy’ was not declared in this scope
 manripper.cpp:78:27: error: ‘strlen’ was not declared in this scope
 manripper.cpp:81:25: error: ‘memset’ was not declared in this scope
 manripper.cpp:82:41: error: ‘memcpy’ was not declared in this scope
 make[2]: *** [manripper.o] Error 1
 make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/manmake'
 make[1]: *** [all-recursive] Error 1
 make[1]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63'
 make: *** [all] Error 2
 kjetil@kjetil:~/yacas/yacas-1.0.63$

 Kjetil




 --
 Statistics  Software Consulting
 GKX Group, GKX Associates Inc.
 tel: 1-877-GKX-GROUP
 email: ggrothendieck at 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.



-- 
Statistics  Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at 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] quantmod getOptionChain Not Work

2012-04-09 Thread R. Michael Weylandt
On a Unix-alike system, if you have an svn client installed, it
suffices to type

svn checkout svn://svn.r-forge.r-project.org/svnroot/quantmod/

at a command prompt -- then use R CMD install on the pkg/ directory to
install the new version. No idea how to do it on Windows

Michael

On Sun, Apr 8, 2012 at 7:11 PM, Sparks, John James jspa...@uic.edu wrote:
 Michael,

 I have not had time to look at this for a while but still wanted to say
 thanks for looking into it and sending this solution.

 By the way, Jeff mentioned that the version of quantmod on the SVN
 (0.3.18) works for this.  I tried to figure out how to download that
 version, but found the documentation on SVN's quite confusing.  Is there
 anyway that you could make that version available?

 Much appreciated.
 --John Sparks



 On Fri, March 23, 2012 5:55 pm, R. Michael Weylandt wrote:
 Sorry about that: two small mistakes and I imagine there are a few
 more I've missed.  This should actually work:

 ###


 library(XML)

 readYahooOptions - function(Symbols, Exp, ...){
   parse.expiry - function(x) {
     if(is.null(x))
       return(NULL)

     if(inherits(x, Date) || inherits(x, POSIXt))
       return(format(x, %Y-%m))

     if (nchar(x) == 5L) {
       x - sprintf(substring(x, 4, 5), match(substring(x,
                                                        1, 3),
 month.abb), fmt = 20%s-%02i)
     }
     else if (nchar(x) == 6L) {
       x - paste(substring(x, 1, 4), substring(x, 5, 6),
                  sep = -)
     }
     return(x)
   }

   clean.opt.table - function(tableIn){
     tableOut - sapply(tableIn[,-2], function(x)
 as.numeric(gsub(,,,x)))
     rownames(tableOut) - tableIn[,2]
     tableOut
   }

   if(missing(Exp))
     optURL -
 paste(paste(http://finance.yahoo.com/q/op?s,Symbols,sep==;),Options,sep=+)
   else
     optURL -
 paste(paste(http://finance.yahoo.com/q/op?s=,Symbols,m=,parse.expiry(Exp),sep=),Options,sep=+)

   if(!missing(Exp)  is.null(Exp)) {
     optPage - readLines(optURL)
     optPage - optPage[grep(View By Expiration, optPage)]
     allExp - gregexpr(m=, optPage)[[1]][-1] + 2
     allExp - substring(optPage, allExp, allExp + 6)
     allExp - allExp[seq_len(length(allExp)-1)] # Last one seems useless ?
     return(structure(lapply(allExp, readYahooOptions,
 Symbols=Symbols), .Names=format(as.yearmon(allExp
   }

   stopifnot(require(XML))

   optURL - readHTMLTable(optURL)

   # Not smart to hard code these but it's a 'good-enough' hack for now
   # Also, what is table 9 on this page?

   list(calls = clean.opt.table(optURL[[10]]),
        puts = clean.opt.table(optURL[[14]]),
        symbol = Symbols)
 }



 On Fri, Mar 23, 2012 at 6:44 PM, R. Michael Weylandt
 michael.weyla...@gmail.com wrote:
 I just got around to taking a look at this, but below is a fix. It
 seems like yahoo finance redesigned the page and rather than reparsing
 all their HTML, I'll use Duncan TL's XML package to make life happier.
 (I loathe HTML parsing)

 This isn't thoroughly tested and it'll break if yahoo redesigns things
 again (I hardcode the table numbers for now) but it seems to work well
 enough. Let me know if you have any errors with it. If Jeff likes it,
 it should be a drop-in replacement for the getOptionChain.yahoo for
 quantmod with a name change.

 Feedback welcome,

 Michael

 #

 library(XML)

 readYahooOptions - function(Symbols, Exp, ...){
 áparse.expiry - function(x) {
 á áif(is.null(x))
 á á áreturn(NULL)

 á áif(inherits(x, Date) || inherits(x, POSIXt))
 á á áreturn(format(x, %Y-%m))

 á áif (nchar(x) == 5L) {
 á á áx - sprintf(substring(x, 4, 5), match(substring(x,
 á á á á á á á á á á á á á á á á á á á á á á á á á á á 1, 3),
 month.abb), fmt = 20%s-%02i)
 á á}
 á áelse if (nchar(x) == 6L) {
 á á áx - paste(substring(x, 1, 4), substring(x, 5, 6),
 á á á á á á á á sep = -)
 á á}
 á áreturn(x)
 á}

 áclean.opt.table - function(tableIn){
 á átableOut - lapply(tableIn[,-2], function(x)
 as.numeric(gsub(,,,x)))
 á árownames(tableOut) - tableIn[,2]
 á}

 áif(missing(Exp))
 á áoptURL -
 paste(paste(http://finance.yahoo.com/q/op?s,Symbols,sep==;),Options,sep=+)
 áelse
 á áoptURL -
 paste(paste(http://finance.yahoo.com/q/op?s=,Symbols,m=,parse.expiry(Exp),sep=),Options,sep=+)

 áif(!missing(Exp)  is.null(Exp)) {
 á áoptPage - readLines(optURL)
 á áoptPage - optPage[grep(View By Expiration, optPage)]
 á áallExp - gregexpr(m=, optPage)[[1]][-1] + 2
 á áallExp - substring(optPage, allExp, allExp + 6)
 á áallExp - allExp[seq_len(length(allExp)-1)] # Last one seems
 useless ? Always true?
 á áreturn(structure(lapply(allExp, readYahooOptions,
 Symbols=Symbols), .Names=format(as.yearmon(allExp
 á}

 ástopifnot(require(XML))

 áoptURL - readHTMLTable(optURL)

 á# Not smart to hard code these but it's a 'good-enough' hack for now
 á# Also, what is table 9 on this page?
 áCALLS - optURL[[10]]
 áPUTS - optURL[[14]]

 álist(calls = 

Re: [R] how to convert seconds to 12 hour time format

2012-04-09 Thread Henrique Dallazuanna
Try this:

Sys.setlocale(category = LC_TIME, locale = US)
format(as.POSIXct('0001-01-01 00:00:00') + x, %I:%M:%S %p)


On Mon, Apr 9, 2012 at 9:38 PM, cassie jones cassiejone...@gmail.com wrote:

 Hello everyone,

 I am wondering if there is any routine in R which can convert time given in
 'seconds' unit to the 12 hour time format. For example, suppose the data
 set looks like

 x=c(36885,84000,20)  #x in seconds

 I want to get the output as

 [1]  11:14:45 AM
 [2]  11:20:00 PM
 [3] 12:20:00 AM

 Does anyone have any idea? Thanks in advance.


 Cassie

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




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

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


Re: [R] how to convert seconds to 12 hour time format

2012-04-09 Thread cassie jones
Thanks Henrique.It worked!!! I appreciate your help.

Cassie

On Mon, Apr 9, 2012 at 8:06 PM, Henrique Dallazuanna www...@gmail.comwrote:

 Try this:

 Sys.setlocale(category = LC_TIME, locale = US)
 format(as.POSIXct('0001-01-01 00:00:00') + x, %I:%M:%S %p)


 On Mon, Apr 9, 2012 at 9:38 PM, cassie jones cassiejone...@gmail.com
 wrote:
 
  Hello everyone,
 
  I am wondering if there is any routine in R which can convert time given
 in
  'seconds' unit to the 12 hour time format. For example, suppose the data
  set looks like
 
  x=c(36885,84000,20)  #x in seconds
 
  I want to get the output as
 
  [1]  11:14:45 AM
  [2]  11:20:00 PM
  [3] 12:20:00 AM
 
  Does anyone have any idea? Thanks in advance.
 
 
  Cassie
 
 [[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.




 --
 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] Problems with CRAN package Ryacas

2012-04-09 Thread Kjetil Halvorsen
see below!

On Mon, Apr 9, 2012 at 7:38 PM, Gabor Grothendieck
ggrothendi...@gmail.com wrote:
 We have had Linux users that successfully used Ryacas but if you can't
 get yacas to build on your particular system then you could try
 running the Windows version under wine.  A Windows binary of yacas is
 available so you won't have to build yacas.  I don't know if anyone
 has tried that yet but its worth a try.

 Another possibility is to try rsympy or rmathpiper (see introductory
 paragraphs on Ryacas home page for links).

 On Mon, Apr 9, 2012 at 8:00 PM, Kjetil Halvorsen
 kjetilbrinchmannhalvor...@gmail.com wrote:
 see below!

 On Mon, Apr 9, 2012 at 3:04 PM, Gabor Grothendieck
 ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 4:02 PM, Gabor Grothendieck
 ggrothendi...@gmail.com wrote:
 On Mon, Apr 9, 2012 at 2:47 PM, Martin Maechler
 maech...@stat.math.ethz.ch wrote:
 Apropos:

 I don't have the problems, the OP had, but on my ubuntu
 notebook, Ryacas does not return expressions (just the strings),
 and hence

    as.expression( yacas-result )

 always gives NULL  and e.g. the   demo(Ryacas-Function)
 also fails:

   yacas(expression(deriv(BurrCDF(x,c,k
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   yy - yacas(expression(deriv(BurrCDF(x,c,k
   yy
  k*c*x^(c-1)*(x^c+1)^(-(k+1));
   str(yy)
  List of 2
  $          : NULL
  $ YacasForm: chr k*c*x^(c-1)*(x^c+1)^(-(k+1));
  - attr(*, class)= chr yacas
  

 I have

 - yacas 1.2.2  (standard ubuntu package)

 The latest version of Ryacas on CRAN is 0.2-11
 http://cran.r-project.org/package=Ryacas

 and the version of yacas supported by Ryacas is:

 yacas(Version())
 expression(1.0.63)

 The troubleshooting section on the home page that I referred to does
 mention this.

I am trying out rsympy by now, and that is working just fine!

Thanks.

Kjetil


 Thanks!

 I deinstalled the synaptic-installed yacas, downloaded yacas_1.0.63.tgz
 unpacked it , run
 ./configure --enable-server
 then did make

 But make does not terminate, it ends in error! I give the last few
 lines of output:

 rm -fr .libs/filescanner.la .libs/filescanner.* .libs/filescanner.*
 gcc -shared  filescanner.lo plugin.lo   -Wl,-soname -Wl,filescanner.so
 -o .libs/filescanner.so
 ar cru .libs/filescanner.a  filescanner.o plugin.o
 ranlib .libs/filescanner.a
 creating filescanner.la
 (cd .libs  rm -f filescanner.la  ln -s ../filescanner.la filescanner.la)
 make[3]: Leaving directory 
 `/home/kjetil/yacas/yacas-1.0.63/plugins/filescanner'
 make[3]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
 make[3]: Nothing to be done for `all-am'.
 make[3]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
 make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/plugins'
 Making all in proteus
 make[2]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/proteus'
 make[2]: Nothing to be done for `all'.
 make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/proteus'
 Making all in manmake
 make[2]: Entering directory `/home/kjetil/yacas/yacas-1.0.63/manmake'
 source='manripper.cpp' object='manripper.o' libtool=no \
        depfile='.deps/manripper.Po' tmpdepfile='.deps/manripper.TPo' \
        depmode=gcc3 /bin/bash ../depcomp \
        g++ -DHAVE_CONFIG_H -I. -I. -I..     -g -O2 -Wall  -c -o manripper.o
 `test -f 'manripper.cpp' || echo './'`manripper.cpp
 manripper.cpp: In function ‘void GetBf(char*, int, FILE*)’:
 manripper.cpp:20:21: error: ‘strlen’ was not declared in this scope
 manripper.cpp: In function ‘void ProcessFile(char*)’:
 manripper.cpp:40:23: error: ‘strlen’ was not declared in this scope
 manripper.cpp:43:30: error: ‘strncmp’ was not declared in this scope
 manripper.cpp:57:25: error: ‘memset’ was not declared in this scope
 manripper.cpp:58:39: error: ‘memcpy’ was not declared in this scope
 manripper.cpp:78:27: error: ‘strlen’ was not declared in this scope
 manripper.cpp:81:25: error: ‘memset’ was not declared in this scope
 manripper.cpp:82:41: error: ‘memcpy’ was not declared in this scope
 make[2]: *** [manripper.o] Error 1
 make[2]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63/manmake'
 make[1]: *** [all-recursive] Error 1
 make[1]: Leaving directory `/home/kjetil/yacas/yacas-1.0.63'
 make: *** [all] Error 2
 kjetil@kjetil:~/yacas/yacas-1.0.63$

 Kjetil




 --
 Statistics  Software Consulting
 GKX Group, GKX Associates Inc.
 tel: 1-877-GKX-GROUP
 email: ggrothendieck at 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.



 --
 Statistics  Software Consulting
 GKX Group, GKX Associates Inc.
 tel: 1-877-GKX-GROUP
 email: ggrothendieck at gmail.com

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

Re: [R] binned tabulation

2012-04-09 Thread Rui Barradas
Hello,


Jean V Adams wrote
 
 Delia,
 
 name - data.frame(Behavior=c(1, 2, 1, 2, 1), Time=c(0, 40, 45, 55, 57))
 
 appear - name$Time[name$Behavior==1]
 disappear - name$Time[name$Behavior==2]
 if(length(appear)  length(disappear)) disappear - c(disappear, 60)
 sum(disappear - appear)
 
 Delia Shelton wrote on 04/09/2012 12:30:23 PM:
 
 Hi,
 
 I am attempting to tabulate binned data. The '1' represents the 
 appearance of the focal mouse pup, and '2' represents the 
 disappearance of the focal mouse pup. The code written below is 
 intended to calculate the total time spent appeared out of 3600s. 
 For Sample 1, both the hand calculation and R code yield the same 
 result, 50. A problem seems to occur when '1' is the last entry. For
 Sample 2, the total time appeared is 53  (hand calculation), 
 however, using the R code below yields 55. If you have any 
 suggestions for solving the problem, please let me know. 
 
 Thank you in advance for any assistance you may provide.
 
 Delia
 
 
 Sample 1
 
 0.0 1
 40 2
 45 1
 55 2
 
 
 Sample 2
 
 0.0 1
 40 2
 45 1
 55 2
 57 1
 
 
 
 name = read.table(file.choose(),header=F) # opening a data file
 colnames(name)-c(Time, Behavior)
 name = data.frame(name$Behavior, name$Time)
 colnames(name)-c(Behavior, Time)
 name-name[name$Time  3600, ];
 
 ## run file 
 
 x-seq(0,3600, by = 60) # total time partition by time which is 60
 
 if (tail(name$Behavior, 1) == 1) {name-rbind(name, c(4, 3600))} else
 {name-rbind(name, c(1, 3600))} 
 
 if (nrow(name) %% 2 != 0)
  {name -name[-c(nrow(name)), -c(nrow(name))]}
 
 q-NULL
 for (y in (1: nrow(name)))
 {
 if (y %% 2 != 0)
 q-c(q, (c(name$Time[y]:name$Time[y +1])))
 }
 
 b-table(cut(q,x))
 
 sum(b)
 
   [[alternative HTML version deleted]]
 
 __
 R-help@ mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide
 http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.
 


Or use a logical index into the data.frame's rows.


fun.count - function(x, increment=60){
if (x$Behavior[nrow(x)] == 1)
x - rbind(x, c(4, increment*ceiling(x$Time[nrow(x)] / 
increment)))
inx - x$Behavior == 1
sum(x$Time[!inx] - x$Time[inx])
}

fun.count(name)

Hope this helps,

Rui Barradas


--
View this message in context: 
http://r.789695.n4.nabble.com/binned-tabulation-tp4543405p4544257.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] Panel.abline would not show beyond a certain slope value

2012-04-09 Thread Aziz, Muhammad Fayez

Yup. Latest version 2.15.0 for windows solved the problem alright! Thank you 
Uwe.

Regards,
Fayez

From: Uwe Ligges [lig...@statistik.tu-dortmund.de]
Sent: Monday, April 09, 2012 9:56 AM
To: Aziz, Muhammad Fayez
Cc: r-help@r-project.org
Subject: Re: [R] Panel.abline would not show beyond a certain slope value

On 09.04.2012 06:31, Aziz, Muhammad Fayez wrote:

 Hi,

 Please find the code and data following. Problem appears in lines 37 - 42 of 
 the code. I am using R.2.13.0 on WinXP.

Works for me with a recent version of R (i.e. 2.15.0) / grid and lattice
(i.e. 0.20-6).
Please always try with recent versions of R and packages before asking.

Best,
Uwe Ligges



 Regards,
 Fayez


 PowerLawGraphsAblineQn.R‎

 library('reshape')
 library(lattice)
 library(igraph) # to use power.law.fit function
 library(latticeExtra) # to use panel.lmlineq in loglog xyplot

 File2Open = C:\\Documents and Settings\\All 
 Users\\Documents\\Desktop\\Fayez\\RPractice\\PowerLawGraphsAblineQnData.txt
 DataTable = read.table(File2Open, header = TRUE, sep = \t)
 md- melt(DataTable, id.vars = c('Range', 'TheValue', 'TheNo')) # removed 
 'DegType' column 030312
 md$variable- as.numeric(substr(md$variable, 9, 
 nchar(as.character(md$variable


 mypanel4loglog-
  function(x, # x is the variable column in melted data, equals the Domain 
 Nos
  y, # y is the value column in melted data, the degrees
  ... # Rest of the arguments
  )
 {

  kfreq- table(y); # compute frquency hash table of y, the values
  k- 1:max(y)
  for (i in k) {
  ichar- as.character(i) # convert to match the names(freq), the 
 character-based hash key of freq, which is the value
  if (!(ichar %in% names(kfreq)))
  kfreq[ichar]- 0
  }
  sortedkeys- as.character(k)
  kfreq-  kfreq[sortedkeys]
  pk- kfreq / length(y)
  logk- log(k)
  logpk- log(pk)

  logpk[logpk == -Inf] =  # remove the -Inf or log(p(k)) = 0 values for 
 lm function, NULL is 0-length so use  instead that has length of one null 
 character
  logpk- as.numeric(logpk) #  converts all values to character, lm 
 needs numeric
  print(rbind(logk, logpk))#write.table(rbind(k, kfreq, pk, logk, 
 logpk), paste(FilePath, \\data, sep=), sep = \t, append = TRUE)

  panel.xyplot(col=blue, logk, logpk, type = c('p', 'r'))

  lm.r = lm(logpk ~ logk)
  panel.abline(coef=c(-4.847634, -1.037480)) # ---  This gets drawn
  panel.abline(coef=c(-4.847634, -1.037481)) # ---  This doesn't get 
 drawn
  print(coef(lm.r)) # -4.847634   -1.349699 ;  -3.377894   -1.498693

 } # end mypanel4loglog

 myprepanel4loglog-
  function(x, # x is the variable column in melted data, equals the 
 Network's ages
  y, # y is the value column in melted data, the degrees
  ... # Rest of the arguments
  )
 {
  FreqTable- as.data.frame(table(y))
  FreqsVector- sort(FreqTable$Freq)
  Min- FreqsVector[1] # first element - the lowest value frequency
 #print(c(Max2, length(y)))
  list(ylim = c(log(Min / length(y)), 0), xlim = c(0, log(max(y # 
 log(p(k)) is always -ve as p(k) is decimal, so max(log(p(k)) is 0
 } # end myprepanel4loglog

 print(xyplot(value ~ variable | Range,
  data = md,
  xlab = log(k); Panel = Range,
  ylab = log(p(k)),
  main = log(k) vs. log(p(k)),
  groups = Range,
  pch = 20, # dots instead of circles
  panel = mypanel4loglog,
  prepanel = myprepanel4loglog, # to set the scale of k and pk
  scales = list(x = (relation = free), y = (relation = free)), # 
 necessary to make x-axis in each panel adjustable according to k
 ))



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
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 Sweave users collaborate with Word users?

2012-04-09 Thread Frank Harrell
See http://biostat.mc.vanderbilt.edu/SweaveConvert for several other
approaches.  The most solid way I've found is to convert PDF to Word.

Frank

Alexander Shenkin wrote
 
 Thanks for the heads up, Paul.  That's good to know.  I happen to be in
 academics, as are my collaborators, but I could imagine running into
 problems down the road.  The good thing seems to be that it's just the
 users who want to interact with R who need the software.  If
 collaborators are just touching the text, or are inserting graphics of
 their own, then sounds like things are still good.  But, that's just my
 impression from the brief description I've seen online.
 
 Thanks,
 Allie
 
 On 4/9/2012 7:42 AM, Paul Bivand wrote:
 If you're considering SWord, you should remember that the licence is
 not the normal R licence and in commercial use will require a
 commercial licence. While some academic disciplines use Word etc, the
 issue may be more common outside academia.

 For those of us where such requirements involve a procurement process,
 the need to purchase something when other users (and the
 administration) are happy with their Word/Excel solutions may be an
 insurmountable barrier.

 Good luck

 Paul Bivand
 Centre for Economic and Social Inclusion (a non-profit organisation)
 London



 On 9 April 2012 13:23, Richard M. Heiberger lt;rmh@gt; wrote:
 You might want to consider SWord, which provides similar facilities for
 the
 Word and R
 user.  Word-oriented co-authors can modify the Word part of the document
 without
 impacting the R part of the document.

 SWord is by Thomas Baier thomas@, author of the statconnDCOM
 interface
 that is underneath RExcel.  See rcom.univie.ac.at for information and
 download and to
 sign up on the rcom email list.

 Rich

 On Sat, Apr 7, 2012 at 4:54 PM, Alexander Shenkin lt;ashenkin@gt;
 wrote:

 Hello All,

 I'm getting my workflow switched over to Sweave, which is very cool.
 However, I collaborate with folks (as many of you must as well) who use
 Word to Track Changes amongst a group while crafting a paper.  In the
 simplest case, there will just be two people (one Sweave user and one
 Word user) editing a paper.

 I'm wondering, how do Sweave users go about this?  I could convert a
 sweave file to a .docx easily enough via an intermediary pdf, rtf, html
 or otherwise.  However, once the file has been marked up with changes,
 the challenge is to migrate those (accepted) changes back to the sweave
 document.  Perhaps the most straightforward way is to manually
 back-propagate changes, but I imagine that could be a painstaking
 process.

 Ideally, I imagine a tool that puts invisible tags in the word document
 when it is originally produced from Sweave, and is then able to
 propagate changes back to that sweave file after markup.  I'd be
 pleasantly surprised if such a tool existed.

 Perhaps there are other ways of making this work.  Any thoughts are
 kindly appreciated!

 Thanks,
 Allie

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

[[alternative HTML version deleted]]

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


-
Frank Harrell
Department of Biostatistics, Vanderbilt University
--
View this message in context: 
http://r.789695.n4.nabble.com/How-do-Sweave-users-collaborate-with-Word-users-tp4540061p4544413.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 add 3d-points to bplot {rms} figure?

2012-04-09 Thread David Winsemius


On Apr 9, 2012, at 6:24 PM, phi771 wrote:


Hello!
I have created a bplot-figure using this code:

*file - 2dcali_red.ttt
ux-as.matrix(read.table(file, dec = ,))

mode(ux)-'numeric'

vel-ux[,1]
ang-ux[,2]
x-ux[,3]
y-ux[,4]

dat- data.frame(ang=ang, x=x,y=y)

require(rms)

ddist2 - datadist(dat)
options(datadist=ddist2)

fitn - lrm(ang ~ rcs(x,4) + rcs(y,4), data=dat)
predi - Predict(fitn, x, y)

bplot(predi,lfun=wireframe, screen = list(z = -40, x = -80),  
drape=TRUE)*


The file 2dcali_red.ttt  consists of 4 columns can be found here :
http://www.color-space.de/upload/dl/2dcali_red.ttt link

The code gives me the bplot-figure, which looks like this:
http://r.789695.n4.nabble.com/file/n4544050/Screen_shot_2012-04-10_at_12.04.24_AM.png

Now I want to add some 3d data (x,y,z) to this plot (like  
scatterplot3d),
but i don't know how. I also don't know how to extract the x,y,z  
mesh as
coordinates from the bplot. Some help would be very appreciated.  
Thank you!


That is most probably a lattice::wireframe object. You could use  
'trellis.focus' or the 'layer' functions in pkg::latticeExtra to  
attempt panel.3dbars() additions. I haven't actually done this, but I  
have added panel function results to the lattice::contourplots that  
are created byrms/Hmisc.


The 3d aspect adds a significant level of extra work: 
http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3274.html

Maybe a contourplot would suffice?

--
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] How to get the SS and MS from oneway.test?

2012-04-09 Thread lulumecindy
Hello everyone:
I'm a new member of this group.
I have a question about oneway.test.
When I use anova(lm()) to analysis the
ANOVA,
I can get the information about Sum Sq  and Mean
Sq.
(The R code and the results are as follows.)
 
  anova(lm(BackCalac~factor(Assay),data=Control))
Analysis of Variance Table
Response: BackCalac
  Df  Sum Sq Mean Sq F valuePr(F)
factor(Assay)  4 270.846  67.711  56.219 1.345e-10 ***
Residuals 20  24.088   1.204 

But it default the variances are the same.
If the variances aren't equal. I need to use the
oneway.test method.
Because oneway.test has the option about
var.equal=F.
Here, I have a question about oneway.test,
How can I get SS, and MS information from
oneway.test?
My R code and the results are as follows. 
Thank you very much.  :)

 oneway.test(BackCalac~factor(Assay), var.equal=T,data=Control)
One-way analysis of means
data:  BackCalac and factor(Assay) 
F = 56.2191, num df = 4, denom df = 20, p-value = 1.345e-10
 oneway.test(BackCalac~factor(Assay), var.equal=F,data=Control)
One-way analysis of means (not assuming equal variances)
data:  BackCalac and factor(Assay) 
F = 92.8834, num df = 4.000, denom df = 9.625, p-value = 1.165e-07


--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-get-the-SS-and-MS-from-oneway-test-tp4544417p4544417.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 add 3d-points to bplot {rms} figure?

2012-04-09 Thread David Winsemius


On Apr 9, 2012, at 10:39 PM, David Winsemius wrote:



On Apr 9, 2012, at 6:24 PM, phi771 wrote:


Hello!
I have created a bplot-figure using this code:

*file - 2dcali_red.ttt
ux-as.matrix(read.table(file, dec = ,))

mode(ux)-'numeric'

vel-ux[,1]
ang-ux[,2]
x-ux[,3]
y-ux[,4]

dat- data.frame(ang=ang, x=x,y=y)

require(rms)

ddist2 - datadist(dat)
options(datadist=ddist2)

fitn - lrm(ang ~ rcs(x,4) + rcs(y,4), data=dat)
predi - Predict(fitn, x, y)

bplot(predi,lfun=wireframe, screen = list(z = -40, x = -80),  
drape=TRUE)*


The file 2dcali_red.ttt  consists of 4 columns can be found here :
http://www.color-space.de/upload/dl/2dcali_red.ttt link

The code gives me the bplot-figure, which looks like this:
http://r.789695.n4.nabble.com/file/n4544050/Screen_shot_2012-04-10_at_12.04.24_AM.png

Now I want to add some 3d data (x,y,z) to this plot (like  
scatterplot3d),
but i don't know how. I also don't know how to extract the x,y,z  
mesh as
coordinates from the bplot. Some help would be very appreciated.  
Thank you!


That is most probably a lattice::wireframe object. You could use  
'trellis.focus' or the 'layer' functions in pkg::latticeExtra to  
attempt panel.3dbars() additions. I haven't actually done this, but  
I have added panel function results to the lattice::contourplots  
that are created byrms/Hmisc.


The 3d aspect adds a significant level of extra work: 
http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3274.html


There is also a panel.3d.identify and panel.identify.cloud that might  
be needed:


https://r-forge.r-project.org/scm/viewvc.php/pkg/R/interaction.R?view=markuprevision=675root=lattice

My efforts at adding bars to that plot have failed so far. This does  
not do anything useful:


 trellis.focus(panel, 1, 1)
 panel.3dbars(x=-1.5, y=-.1, z=0,distance=20,  col=red,
 xlim=c(-2.4,-1), ylim=c(-0.15, 0.15), zlim=c(-50,100),
  xlim.scaled=c(-0.5,0.5), ylim.scaled=c(-0.5,0.5),
  zlim.scaled=c(-0.5,0.5), zero.scaled=25)
 trellis.unfocus()



Maybe a contourplot would suffice?




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] how to add 3d-points to bplot {rms} figure?

2012-04-09 Thread Bert Gunter
Below.

-- Bert

On Mon, Apr 9, 2012 at 7:39 PM, David Winsemius dwinsem...@comcast.net wrote:

 On Apr 9, 2012, at 6:24 PM, phi771 wrote:

 Hello!
 I have created a bplot-figure using this code:

 *file - 2dcali_red.ttt
 ux-as.matrix(read.table(file, dec = ,))

 mode(ux)-'numeric'

 vel-ux[,1]
 ang-ux[,2]
 x-ux[,3]
 y-ux[,4]

 dat- data.frame(ang=ang, x=x,y=y)

 require(rms)

 ddist2 - datadist(dat)
 options(datadist=ddist2)

 fitn - lrm(ang ~ rcs(x,4) + rcs(y,4), data=dat)
 predi - Predict(fitn, x, y)

 bplot(predi,lfun=wireframe, screen = list(z = -40, x = -80), drape=TRUE)*

 The file 2dcali_red.ttt  consists of 4 columns can be found here :
 http://www.color-space.de/upload/dl/2dcali_red.ttt link

 The code gives me the bplot-figure, which looks like this:

 http://r.789695.n4.nabble.com/file/n4544050/Screen_shot_2012-04-10_at_12.04.24_AM.png

 Now I want to add some 3d data (x,y,z) to this plot (like scatterplot3d),
 but i don't know how. I also don't know how to extract the x,y,z mesh as
 coordinates from the bplot. Some help would be very appreciated. Thank
 you!


 That is most probably a lattice::wireframe object. You could use
 'trellis.focus' or the 'layer' functions in pkg::latticeExtra to attempt
 panel.3dbars() additions. I haven't actually done this, but I have added
 panel function results to the lattice::contourplots that are created
 byrms/Hmisc.

 The 3d aspect adds a significant level of extra work:
And is generally a bad idea. Perspective plots of all types are much
over-hyped. They are pretty, but the very act of creating artificial
perspective for non-realistic shapes can distort perception and lead
to misinterpretation of data. Parts of the scene hide other parts; it
is difficult to quantitatively judge depths of valleys and heights of
hills; etc. And, as David notes, it is very difficult to get right.

 http://tolstoy.newcastle.edu.au/R/e2/help/06/10/3274.html

 Maybe a contourplot would suffice?
-- and is generally a much better, albeit less sexy, idea.


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



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm

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


Re: [R] How to get the SS and MS from oneway.test?

2012-04-09 Thread Rolf Turner


It is not clear to me how meaningful sums of squares and mean squares
are in the context of oneway.test.().  Be that as it may, the short answer
is you ***can't*** get at them from the output of oneway.test().  The
oneway.test() function simply does not return enough information.

If you ***really*** want/need (or think you do!) SS and MS information,
you'll simply have to hack the code of oneway.test().  The code looks
simple enough.  I would advise reading the cited 1951 Biometrics paper
by Welch in conjunction with your code-hacking.

Good luck.

cheers,

Rolf Turner

On 10/04/12 14:16, lulumecindy wrote:

Hello everyone:
 I'm a new member of this group.
 I have a question about oneway.test.
 When I use anova(lm()) to analysis the
ANOVA,
 I can get the information about Sum Sq  and Mean
Sq.
 (The R code and the results are as follows.)

anova(lm(BackCalac~factor(Assay),data=Control))
Analysis of Variance Table
Response: BackCalac
   Df  Sum Sq Mean Sq F valuePr(F)
factor(Assay)  4 270.846  67.711  56.219 1.345e-10 ***
Residuals 20  24.088   1.204

 But it default the variances are the same.
 If the variances aren't equal. I need to use the
oneway.test method.
 Because oneway.test has the option about
var.equal=F.
 Here, I have a question about oneway.test,
 How can I get SS, and MS information from
oneway.test?
 My R code and the results are as follows.
 Thank you very much.  :)


oneway.test(BackCalac~factor(Assay), var.equal=T,data=Control)

 One-way analysis of means
data:  BackCalac and factor(Assay)
F = 56.2191, num df = 4, denom df = 20, p-value = 1.345e-10

oneway.test(BackCalac~factor(Assay), var.equal=F,data=Control)

 One-way analysis of means (not assuming equal variances)
data:  BackCalac and factor(Assay)
F = 92.8834, num df = 4.000, denom df = 9.625, p-value = 1.165e-07


--
View this message in context: 
http://r.789695.n4.nabble.com/How-to-get-the-SS-and-MS-from-oneway-test-tp4544417p4544417.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] How to get the SS and MS from oneway.test?

2012-04-09 Thread Tal Galili
From what I can see the function does not give back this information
(either in the print output or the output object structure itself).

The thing is that I have never studied the statistic used for this test, so
I am not sure how it works and if giving the SS MS is meaningful or not.
But if you want to dig in, simply write in R the command:
oneway.test

The code that seems to be relevant is near the end:
}
else {
m - sum(w.i * m.i)/sum.w.i
STATISTIC - sum(w.i * (m.i - m)^2)/((k - 1) * (1 + 2 *
(k - 2) * tmp))
PARAMETER - c(k - 1, 1/(3 * tmp))
PVAL - pf(STATISTIC, k - 1, 1/(3 * tmp), lower.tail = FALSE)
METHOD - paste(METHOD, (not assuming equal variances))
}

So you could simply copy paste the full code of the function and edit it to
extract this information...







Contact
Details:---
Contact me: tal.gal...@gmail.com |  972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
www.r-statistics.com (English)
--




On Tue, Apr 10, 2012 at 5:16 AM, lulumecindy cindy852...@gmail.com wrote:

 Hello everyone:
I'm a new member of this group.
I have a question about oneway.test.
When I use anova(lm()) to analysis the
 ANOVA,
I can get the information about Sum Sq  and Mean
 Sq.
(The R code and the results are as follows.)

   anova(lm(BackCalac~factor(Assay),data=Control))
 Analysis of Variance Table
 Response: BackCalac
  Df  Sum Sq Mean Sq F valuePr(F)
 factor(Assay)  4 270.846  67.711  56.219 1.345e-10 ***
 Residuals 20  24.088   1.204

But it default the variances are the same.
If the variances aren't equal. I need to use the
 oneway.test method.
Because oneway.test has the option about
 var.equal=F.
Here, I have a question about oneway.test,
How can I get SS, and MS information from
 oneway.test?
My R code and the results are as follows.
Thank you very much.  :)

  oneway.test(BackCalac~factor(Assay), var.equal=T,data=Control)
One-way analysis of means
 data:  BackCalac and factor(Assay)
 F = 56.2191, num df = 4, denom df = 20, p-value = 1.345e-10
  oneway.test(BackCalac~factor(Assay), var.equal=F,data=Control)
One-way analysis of means (not assuming equal variances)
 data:  BackCalac and factor(Assay)
 F = 92.8834, num df = 4.000, denom df = 9.625, p-value = 1.165e-07


 --
 View this message in context:
 http://r.789695.n4.nabble.com/How-to-get-the-SS-and-MS-from-oneway-test-tp4544417p4544417.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.