Re: [R] [FORGED] How to alpha entire plot?

2018-06-04 Thread Paul Murrell

Hi

Here is one way to do it ...

col1 <- "red"
col2 <- "blue"
EU <- data.frame(EuStockMarkets)
with(EU, plot(DAX, CAC, col=col2, type="h", ylim=c(0,6000)))
par(new=TRUE)
with(EU, plot(DAX, FTSE, col=col1, type="h", ylim=c(0,6000)))

## Convert 'graphics' to 'grid'
library(gridGraphics)
grid.echo()
## grid.ls(print=grobPathListing, viewports=TRUE)

## Rasterize spikes
library(rasterize)
downViewport("graphics-window-1-1")
grid.rasterize("graphics-plot-1-spike-1")
upViewport(0)
downViewport("graphics-window-2-1")
grid.rasterize("graphics-plot-2-spike-1")
upViewport(0)

## Apply alpha adjustment to rasterized spikes
spike1 <- as.matrix(grid.get("graphics-plot-1-spike-1")$raster)
alphaSpike1 <- adjustcolor(spike1, alpha=.3)
dim(alphaSpike1) <- dim(spike1)
grid.edit("graphics-plot-1-spike-1", raster=as.raster(alphaSpike1))
spike2 <- as.matrix(grid.get("graphics-plot-2-spike-1")$raster)
alphaSpike2 <- adjustcolor(spike2, alpha=.3)
dim(alphaSpike2) <- dim(spike2)
grid.edit("graphics-plot-2-spike-1", raster=as.raster(alphaSpike2))

... though that may not be the best way and may just reflect what I have 
been thinking about recently ...


https://www.stat.auckland.ac.nz/~paul/Reports/rasterize/rasterize.html

Paul


On 01/06/18 08:56, Ed Siefker wrote:

I have two chromatograms I want plotted on the same axes.
I would like the plots to be transparent, so the first chart is
not obscured.

I have tried adjustcolor(..., alpha.f=0.3), the problem is that
my chromatogram is so dense with datapoints that they
overlap and the entire graph just ends up a solid color.  The
second histogram still obscures the first.

Consider this example:


col1 <- adjustcolor("red", alpha.f=0.3)
col2 <- adjustcolor("blue", alpha.f=0.3)
EU <- data.frame(EuStockMarkets)
with(EU, plot(DAX, CAC, col=col2, type="h", ylim=c(0,6000)))
par(new=TRUE)
with(EU, plot(DAX, FTSE, col=col1, type="h", ylim=c(0,6000)))

The density of the red plot around 2000 completely obscures the blue
plot behind it.

What I would like to do is plot both plots in solid colors, then alpha
the entire thing, and then overlay them.  Or some other method that
achieves a comparable result.
Thanks

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



--
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag 92019
Auckland
New Zealand
64 9 3737599 x85392
p...@stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/

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


Re: [R] Running segmented on grouped data and collating model parameters in a data frame

2018-06-04 Thread Vito M. R. Muggeo

dear Sumi,
I am not familiar with the dplyr package (%>%..), however if you want to 
fit the model for each subject times freq interaction, a simple for loop 
will suffice.

Here possible code:

Assuming d is the dataframe, something like

subj<-levels(d$subject)
fr<-unique(d$freq)

#new dataframe where the estimates will be stored
newd<-expand.grid(subj=subj, freq=fr)
newd<-cbind(newd, matrix(-99,nrow(newd),9))
names(newd)[-(1:2)]<-c(paste0("slope",1:3), 
paste0(c("CI.inf.","CI.sup."), rep(paste0("slope",1:3),each=2)))


library(segmented)

#fit the model for each subject x freq combination
for(i in subj){
  for(j in fr){
current.d<-subset(d, subject==i & freq==j)
if(nrow(current.d)>0){
  o<-lm(Ldp_p~L2, data=current.d)
  os<-try(segmented(o, ~L2, psi=c(20,50)), silent=TRUE)
  if(class(os)[1]!="try-error"){
 est.slope<-as.vector(slope(os)[[1]][,1]) #point est
 ci.slope<-as.vector(t(slope(os)[[1]][,4:5])) #CI
 newd[newd$subj==i & newd$freq==j, 
-c(1:2)]<-c(est.slope, ci.slope)

   }
 }

  }
}


However it seems that you can be interested in segmented *mixed* 
models.. I mean, rather than estimating a segmented model for each 
subject, you could estimate a single model assuming random effects (for 
each model parameter, including the breakpoint) for each subject:


Have a look to

http://journals.sagepub.com/doi/abs/10.1177/1471082X13504721

Let me know (not on the R list) if you are interested in relevant code

best,
vito


Il 03/06/2018 18:03, Sumitrajit Dhar ha scritto:

Hi Folks,

I am trying to teach myself how to solve the problem described below but am
running out of time. Hence the plea for help. Thanks in advance.

Here is my data frame.


t

# A tibble: 12 x 12
subject  ageGrp ear   hearingGrp sexfreqL2   Ldp PhidpNF
SNR   Ldp_p
  
   
  1 HALAF573 A  L A  F 2 0 -19.6 197.  -28.5
8.88  2.10e-6
  2 HALAF573 A  L A  F 2 2 -18.7 203.  -22.0
3.25  2.31e-6
  3 HALAF573 A  L A  F 2 4 -29.1 255.  -27.4
  -1.64  7.04e-7
  4 HALAF573 A  L A  F 2 6 -12.4 174.  -12.2
  -0.206 4.78e-6
  5 HALAF573 A  L A  F 4 0 -28.6 232.  -26.7
  -1.87  7.45e-7
  6 HALAF573 A  L A  F 4 2 -27.2 351.  -28.8
1.59  8.68e-7
  7 HALAF573 A  L A  F 4 4 -20.4  26.2 -35.0
  14.6   1.92e-6
  8 HALAF573 A  L A  F 4 6 -20.0  85.1 -29.8
9.75  2.00e-6
  9 HALAF573 A  L A  F 8 0 -22.8  39.2 -22.1
  -0.689 1.45e-6
10 HALAF573 A  L A  F 8 2 -14.5  13.4 -20.7
6.23  3.76e-6
11 HALAF573 A  L A  F 8 4 -17.3 345.  -21.6
4.30  2.73e-6
12 HALAF573 A  L A  F 8 6 -14.1 320.  -21.7
7.59  3.94e-6

# Note that there are more levels of L2 (31 in total)  and 344 other
subjects but I truncated the frame for posting.

# I want to do this:  t %>%  group_by(freq) %>% [run segmented] %>% [create
a data frame with [subject, freq, breakpoint1, breakpoint2, slope1, slope2,
slope3, L2 when Ldp_p == 0].

# Also note that ultimately I will be grouping by "subject, freq".

# I can run the models and get believable results. The following run on a
data frame with L2 between 0 and 60.

out.lm <- lm(Ldp_p ~ L2, data = t)
o <- segmented(out.lm, seg.Z = ~L2, psi = list(L2 = c(20,45)),
control = seg.control(display = FALSE)

o$psi
#Initial Est.St.Err
#psi1.L2  20 30.78256 0.5085192
#psi2.L2  45 53.16390 0.4671701

slope(o)
slope(o)
$L2
#  Est.St.Err. t value   CI(95%).l   CI(95%).u
#slope1  1.2060e-06 1.6606e-07  7.2622  8.6397e-07  1.5480e-06
#slope2  1.0708e-05 2.9196e-07 36.6770  1.0107e-05  1.1309e-05
#slope3 -4.5791e-06 1.3694e-06 -3.3439 -7.3995e-06 -1.7588e-06

Thanks again for bailing me out.

Regards,
Sumit

[[alternative HTML version deleted]]

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



--
==
Vito M.R. Muggeo
Dip.to Sc Econom, Az e Statistiche
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 23895240
fax: 091 485726
http://dssm.unipa.it/vmuggeo
Associate Editor, Statistical Modelling
Chair, Statistical Modelling Society

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

[R-es] rseek 2?

2018-06-04 Thread Javier Marcuzzi
Estimados

Estoy realizando una referencia y veo que rseek antes me daba una cantidad
de resultados a la consulta, hoy la cantidad de resultados es menor, pero
dice en el logo rseek2.

¿Cambió el algoritmo?, ¿es diferente?, ¿se puede acceder al anterior?

¿Alguno sabe que paso? El dato importante para conocer es la cantidad de
resultados que da ante determinada palabra, pero esto me cambio (en números
grandes).

Javier Rubén Marcuzzi

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es


Re: [R] Time and date conversion

2018-06-04 Thread peter dalgaard



> On 4 Jun 2018, at 10:45 , Christofer Bogaso  
> wrote:
> 
> Hi,
> 
> I have an automatic data feed and I obtained a Date vector in the following
> format:
> 
>> Date
> [1] "03 Jun 2018 10:01 am CT""01 Jun 2018 22:04:25 pm CT"
> 
> I now like to convert it to UTC time-zone
> 
> Is there any easy way to convert them so, particularly since 1st element
> doesnt have any Second element whereas the 2nd element has.

..and it also mixes up am/pm notation and 24hr clock.

There are two basic approaches to the format inconsistency thing:

(A) preprocess using gsub() constructions 

> gsub(" (..:..) ", " \\1:00 ", d.txt)
[1] "03 Jun 2018 10:01:00 am CT" "01 Jun 2018 22:04:25 pm CT"

(B) Try multiple formats

> d <- strptime(d.txt, format="%d %B %Y %H:%M:%S %p")
> d[is.na(d)] <- strptime(d.txt[is.na(d)], format="%d %B %Y %H:%M %p")
> d
[1] "2018-06-03 10:01:00 CEST" "2018-06-01 22:04:25 CEST"

I would likely go for (A) since you probably need to do something gsub-ish to 
get the TZ thing in place.

-pd

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

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd@cbs.dk  Priv: pda...@gmail.com

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


Re: [R] Time-series moving average question

2018-06-04 Thread Bill Poling
Good morning Don, thank you for your support. I will give this all more 
investigation and submit my findings when finalized.

Cheers and thanks again for your help Sir.

WHP



From: MacQueen, Don [mailto:macque...@llnl.gov]
Sent: Friday, June 01, 2018 5:03 PM
To: Bill Poling ; r-help (r-help@r-project.org) 

Subject: Re: [R] Time-series moving average question

For your first question, you’re doing moving averages of 3 points (I assume 
that’s what order=3 does). For any given time point of your input data, that 
would be one before, one at, and one after the given time point. Do all of  
your input times have both one before and one after?

Don’t know about your same 8 point forecast values question, not without 
running it myself, but I would say that if the data really is different, yet 
the forecasts are the same, it would have to be coincidental. But a priori this 
seems unlikely, so I’d inspect your script very carefully for mistakes.

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509



From: Bill Poling mailto:bill.pol...@zelis.com>>
Date: Friday, June 1, 2018 at 10:43 AM
To: "MacQueen, Don" mailto:macque...@llnl.gov>>, array 
R-help mailto:r-help@r-project.org>>
Subject: RE: [R] Time-series moving average question

Hi Don, wow, you are so right. I picked that piece up from the bloggers 
tutorial and since I am R naive yet, I thought it was all one step

moving_average = forecast(ma(tdat[1:31], order=2), h=5)

Truly, I usually print and check at every step I can, as painful as it is 
sometimes.
Great lesson for this novice usR.

So the first and last values are NA in each case? Do you know why? Should I 
replace the NA with a value, say the average of the others?

Also, I have 5 series of data I am working with using this script and of course 
each gave me that warning, but only the one series had the same 8 Point 
Forecast values, is that coincidental you think?

Terrific of you to help, I really appreciate it.

WHP


From: MacQueen, Don [mailto:macque...@llnl.gov]
Sent: Friday, June 01, 2018 12:54 PM
To: Bill Poling mailto:bill.pol...@zelis.com>>; r-help 
(r-help@r-project.org) 
mailto:r-help@r-project.org>>
Subject: Re: [R] Time-series moving average question

You are right that there are no NAs in the practice data. But there are NAs in 
the moving average data.

To see this, break your work into two separate steps, like this:

tnr.ma <- ma(dat3[1:28], order=3)
TNR_moving_average <- forecast(tnr.ma, h=8)

I think you will find that the warning comes from the second step.

Print tnr.ma and you will see some NAs.

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509



From: Bill Poling mailto:bill.pol...@zelis.com>>
Date: Friday, June 1, 2018 at 8:58 AM
To: "MacQueen, Don" mailto:macque...@llnl.gov>>, array 
R-help mailto:r-help@r-project.org>>
Subject: RE: [R] Time-series moving average question

Hello Don, thank you for your response. I appreciate your help.

I am using the forecast package, originally I found it following a forecasting 
example on bloggers.com

https://www.r-bloggers.com/time-series-analysis-using-r-forecast-package/

And subsequently located the complete pdf 
https://cran.r-project.org/web/packages/forecast/forecast.pdf

Since I created this practice data using the structure() function I am unsure 
why there would be NA’s as there are none apparently in the structure?

No worries though, I am going to reach out to the package author.

Cheers.

WHP

From: MacQueen, Don [mailto:macque...@llnl.gov]
Sent: Friday, June 01, 2018 11:24 AM
To: Bill Poling mailto:bill.pol...@zelis.com>>; r-help 
(r-help@r-project.org) 
mailto:r-help@r-project.org>>
Subject: Re: [R] Time-series moving average question

My guess would be that if you inspect the output from
ma(dat3[1:28], order=3)
you will find some NAs in it. And then forecast() doesn't like NAs.

But I can't check, because I can't find the ma() and forecast() functions. I 
assume they come from some package you installed; it would be helpful to say 
which package.

-Don

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509


__
mailto:R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] Time and date conversion

2018-06-04 Thread Christofer Bogaso
Hi,

I have an automatic data feed and I obtained a Date vector in the following
format:

> Date
[1] "03 Jun 2018 10:01 am CT""01 Jun 2018 22:04:25 pm CT"

I now like to convert it to UTC time-zone

Is there any easy way to convert them so, particularly since 1st element
doesnt have any Second element whereas the 2nd element has.

Thanks for any pointer.

[[alternative HTML version deleted]]

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


Re: [R] package colorspace and .WhitePoint question

2018-06-04 Thread Achim Zeileis

Glenn,

currently, this is currently not exposed in "colorspace" AFAICS. You can 
modify it by changing .WhitePoint inside colorspace's NAMESPACE, though:


R> assignInNamespace(".WhitePoint", rbind(c(95, 100, 105)),
+ns = "colorspace")
R> as(XYZ(100, 100, 100), "LAB")
   LAB
[1,] 100 8.622384 3.226371

I'll have another look whether this could be exposed easily (cc also 
Paul).


Best,
Z

On Mon, 4 Jun 2018, Glenn Davis wrote:


In colorspace.R  I see:

   setAs("color", "LAB", function(from)
 LAB(.Call("as_LAB", from@coords, class(from), .WhitePoint, PACKAGE = "
colorspace"),
 names = dimnames(from@coords)[[1]]))
   ...
   .WhitePoint = NULL

and then in colorspace.c and the function CheckWhite(),
I see that .WhitePoint = NULL is converted to D65.

I would like to pass a different .WhitePoint to
   as( XYZ( 100,100,100)  , "LAB" )


I have tried 3 methods:
   as( XYZ( 100,100,100)  , "LAB", .WhitePoint=XYZ(95,100,105) )
   .WhitePoint = XYZ(95,100,105)
   assign( ".WhitePoint", XYZ(95,100,105), env=as.environment('package:
colorspace') )
but all fail, for different reasons.

How can I transform XYZ to LAB using a whitepoint different than D65 ?

Thanks,
Glenn Davis
gda...@gluonics.com

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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 -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/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] package colorspace and .WhitePoint question

2018-06-04 Thread Glenn Davis
 In colorspace.R  I see:

setAs("color", "LAB", function(from)
  LAB(.Call("as_LAB", from@coords, class(from), .WhitePoint, PACKAGE = "
colorspace"),
  names = dimnames(from@coords)[[1]]))
...
.WhitePoint = NULL

and then in colorspace.c and the function CheckWhite(),
I see that .WhitePoint = NULL is converted to D65.

I would like to pass a different .WhitePoint to
as( XYZ( 100,100,100)  , "LAB" )


I have tried 3 methods:
as( XYZ( 100,100,100)  , "LAB", .WhitePoint=XYZ(95,100,105) )
.WhitePoint = XYZ(95,100,105)
assign( ".WhitePoint", XYZ(95,100,105), env=as.environment('package:
colorspace') )
but all fail, for different reasons.

How can I transform XYZ to LAB using a whitepoint different than D65 ?

Thanks,
Glenn Davis
gda...@gluonics.com

[[alternative HTML version deleted]]

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