Como desenhar quantil de banda (em R)
Pergunta
Eu tenho um arquivo CSV que contém linhas para cada (Java GC) Evento que me interessa.O objeto consiste em uma subsecond carimbo de data / hora (não eqüidistantes) e algumas variáveis.O objeto se parece com isso:
gcdata <- read.table("http://bernd.eckenfels.net/view/gc1001.ygc.csv",header=TRUE,sep=",", dec=".")
start = as.POSIXct(strptime("2012-01-01 00:00:00", format="%Y-%m-%d %H:%M:%S"))
gcdata.date = gcdata$Timestamp + start
gcdata = gcdata[,2:7] # remove old date col
gcdata=data.frame(date=gcdata.date,gcdata)
str(gcdata)
Resultados
'data.frame': 2997 obs. of 7 variables:
$ date : POSIXct, format: "2012-01-01 00:00:06" "2012-01-01 00:00:06" "2012-01-01 00:00:18" ...
$ Distance.s. : num 0 0.165 11.289 9.029 11.161 ...
$ YGUsedBefore.K.: int 1610619 20140726 20148325 20213304 20310849 20404772 20561918 21115577 21479211 21544930 ...
$ YGUsedAfter.K. : int 7990 15589 80568 178113 272036 429182 982841 1346475 1412181 1355412 ...
$ Promoted.K. : int 0 0 0 0 8226 937 65429 71166 62548 143638 ...
$ YGCapacity.K. : int 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 22649280 ...
$ Pause.s. : num 0.0379 0.022 0.0287 0.0509 0.109 ...
Neste caso, eu me importo com o tempo de Pausa (em segundos).Eu quero desenhar um diagrama, que mostra para cada um (relógio de parede) horas, basicamente, a média como uma linha, a 2% e 98%, como cinza corredor e o valor máximo (a cada hora) como uma linha vermelha.
Eu tenho feito alguns trabalhos, mas usando o q98 funções é feio, ter que usar várias linhas demonstrações parece ser wastefull, e eu não sei como alcançar uma área cinzenta entre o q02 e q98:
q02 <- function(x, ...) { x <- quantile(x,probs=c(0.2)) }
q98 <- function(x, ...) { x <- quantile(x,probs=c(0.98)) }
hours = droplevels(cut(gcdata$date, breaks="hours")) # can I have 2 hours?
plot(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=max),ylim=c(0,2), col="red", ylab="Pause(s)", xlab="Days") # Is always black?
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q98),ylim=c(0,2), col="green")
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q02),ylim=c(0,2), col="green")
lines(aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=mean),ylim=c(0,2), col="blue")
Agora, isso resulta em um gráfico que tem pontos pretos como máximo, uma linha azul, como a média horária e um inferior e superior a 0,2 + 0,98 linha verde.Eu acho que seria melhor legível para ter um corredor cinza, talvez um tracejado máximo (vermelho) linha e de alguma forma, corrigir os rótulos do eixo. Alguma sugestão?(o arquivo está disponível acima)
Solução
bom ver colegas Debian veteranos aqui :) Sua resposta já está muito bom.Como eu trabalho muito com séries de tempo, eu pensei que eu ia jogar uma variante usando o excelente jardim zoológico e xts os pacotes.Este último é construído sobre a antiga e tem, entre outras coisas, a period.apply()
função podemos usar aqui, juntamente com o endpoints()
função para obter a você dois-horária agregados.
Então, no início eu usaria
library(zoo) # for zoo objects
library(xts) # for period.apply
gcdata <- read.table("http://bernd.eckenfels.net/view/gc1001.ygc.csv",
header=TRUE, sep=",", dec=".")
timestamps <- gcdata$Timestamp +
as.POSIXct(strptime("2012-01-01 00:00:00",
format="%Y-%m-%d %H:%M:%S"))
gcdatazoo <- zoo(gcdata[-1], order.by=timestamps) # as zoo object
para criar um zoo
objecto.Sua função permanece:
plotAreaCorridor <- function(x, y, col.poly1="lightgray", col.poly2="gray",...) {
x.pol <- c(x, rev(x), x[1])
y.pol <- c(y[,1], rev(y[,5]),y[,1][1])
plot(x, y[,6]+1, type="n", ...)
polygon(x.pol, y.pol, col=col.poly1, lty=0)
x.pol <- c(x, rev(x), x[1])
y.pol <- c(y[,2], rev(y[,4]), y[,1][1])
polygon(x.pol, y.pol, col=col.poly2, lty=0)
lines(x, y[,3], col="blue") # median
lines(x, y[,6], col="red") # max
invisible(NULL)
}
E em seguida podemos simplificar um pouco:
agg <- period.apply(gcdatazoo[,"Pause.s."], # to which data
INDEX=endpoints(gcdatazoo, "hours", k=2), # every 2 hours
FUN=function(x) quantile(x, # what fun.
probs=c(5,20,50,80,95,100)/100))
#v99 = q99(gcdata$Pause.s.) # what is q99 ?
v99 <- mean(agg[,5]) # mean of 95-th percentile?
plotAreaCorridor(index(agg), # use time index as x axis
coredata(agg), # and matrix part of zoo object as data
ylim=c(0,max(agg[,5])*1.5),
ylab="Quantiles of GC events",
main="NewPar Collection Activity")
abline(h=median(gcdatazoo[,"Pause.s."]), col="lightblue")
abline(h=v99, col="grey")
labeltxt <- paste("99%=",round(v99,digits=3),"s n=", nrow(gcdatazoo),sep="")
text(x=index(agg)[20], y=1.5*v99, labeltxt, col="grey", pos=3) # or legend()
o que dá
O eixo agora é automática e mostra somente nos dias de semana como a amplitude é menor do que na semana;isso pode ser substituir, conforme necessário.
Outras dicas
Você tem que tentar polygon
.Este código pode ser útil:
y98 = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q98)
y02 = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=q02)
ymax = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=max)
ymin = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=min)
ymean = aggregate(gcdata$Pause.s. ~ hours, data=gcdata, FUN=mean)
x = ymean[,1]
y1 = cbind(y02[,2], ymean[,2], y98[,2])
y2 = cbind(ymin[,2], ymean[,2], ymax[,2])
plotAreaCI(x,y2, ylim=c(0,2), xlab="time", ylab="variable")
plotAreaCI(x,y1, ylim=c(0,2), poly.col="blue", add=TRUE)
ou
plotAreaCI(x,y2, ylim=c(0,2), xlab="time", ylab="variable", nice.x = TRUE)
plotAreaCI(x,y1, ylim=c(0,2), mean.lwd=2, poly.col="blue", add=TRUE)
onde a função plotAreaCI
é definido por:
plotAreaCI = function(x, y, add=FALSE, nice.x = FALSE,
xlim=NULL, ylim=NULL,
mean.col="black", mean.lwd=1.5,
poly.col="gray", poly.lty=3,
xlab=NULL, ylab=NULL, main="",
...) {
isFactorX = isClass("factor", x)
if(isFactorX) {
x.label = x
x = as.numeric(x)
}
if(is.null(xlim)) xlim=range(x, na.rm=TRUE)
if(is.null(ylim)) ylim=range(y, na.rm=TRUE)
x.pol = c(x, rev(x), x[1])
y.pol = c(y[,1], rev(y[,3]), y[,1][3])
if(!add) {
plot.new()
plot.window(xlim=xlim, ylim=ylim, ...)
if(!nice.x & isFactorX) {
axis(1, at=x, labels=x.label)
} else {
xticks = axTicks(1)
if(isFactorX) {
xticks = xticks[xticks>=1]
axis(1, at=xticks, labels=x.label[xticks])
} else {
axis(1)
}
}
axis(2, las=1)
box()
title(xlab=xlab, ylab=ylab, main=main)
}
polygon(x.pol, y.pol, col=poly.col, lty=poly.lty)
lines(x, y[,2], col=mean.col, lwd=mean.lwd)
return(invisible())
}
Este é o código que eu uso para desenhar variação temporal do laboratório de analitos (pressão arterial sistólica neste exemplo):
SBP.qtr.mat <- aggregate(set1HLI$SBP,
list( year(set1HLI$Drawdt)+0.25* quarter(set1HLI$Drawdt)),
quantile, prob=c(0.1,0.25,0.5,0.75, 0.9,0.95, 0.975), na.rm=TRUE)
matplot(SBP.qtr.mat[,1], SBP.qtr.mat$x, type="pl")
Não deve ser muito difícil adaptá-lo para o seu problema....ou você poderia postar um exemplo pode ser reproduzido para trabalhar.Com isso, o 10, 25, 50, 75, 90, 95 e 97,5 th percentis em uma única data.quadro e matplot lida com o desenho de um objeto.
A área a cinzento?, ...A abordagem usual é desenhar um polígono indo para o limite inferior, "virada" na extrema direita e chegando de volta no lado de alta, e ligar de volta ao lado esquerdo.O polygon
os argumentos são definidos como x, y
.Há um col
argumento que você deve definir para "cinza".
Para fazer o "2 horas" de sequências para o qual você pode mesclar o dataframe ou utilização comcut.POSIXt" as a breaks argument , there is the option of using multiples of time units with
seq.POSIXt`:
> seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "10 years")
[1] "1910-01-01 12:00:00 GMT" "1920-01-01 12:00:00 GMT" "1930-01-01 12:00:00 GMT" "1940-01-01 12:00:00 GMT"
[5] "1950-01-01 12:00:00 GMT" "1960-01-01 12:00:00 GMT" "1970-01-01 12:00:00 GMT" "1980-01-01 12:00:00 GMT"
[9] "1990-01-01 12:00:00 GMT"
Eu não vi isso documentado, mas você pode usar múltiplos do intervalo com cut.POSIXt
:
> str( cut( seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "years"), "10 years") )
Factor w/ 9 levels "1910-01-01","1920-01-01",..: 1 1 1 1 1 1 1 1 1 1 ...
> str( cut( seq(ISOdate(1910,1,1), ISOdate(1999,1,1), "years"), "5 years") )
Factor w/ 18 levels "1910-01-01","1915-01-01",..: 1 1 1 1 1 2 2 2 2 2 ...
Não estou chegou à seguinte script (ainda precisa de olhar para o mais avançado resposta do DWin).Agora um pouco parece que eu estava procurando, mas o código ainda está muito feio (Por exemplo, eu não sei como alinhar o rótulo e como obter bom xlab rótulos):
plotAreaCorridor = function(x, y, col.poly1="lightgray", col.poly2="gray",...) {
x.pol = c(x, rev(x), x[1])
y.pol = c(y[,1], rev(y[,5]),y[,1][1])
plot(x, y[,6]+1, type="n", ...) # ugly since type="n" does not work for factor
polygon(x.pol, y.pol, col=col.poly1, lty=0)
x.pol = c(x, rev(x), x[1])
y.pol = c(y[,2], rev(y[,4]), y[,1][1])
polygon(x.pol, y.pol, col=col.poly2, lty=0)
lines(x, y[,3], col="blue") # median
lines(x, y[,6], col="red") # max
return(invisible())
}
pause = gcdata$Pause.s.
hours = droplevels(cut(gcdata$date, breaks="hours")) # can I have 2 hours?
agg = aggregate(pause ~ hours, FUN=quantile, probs=c(5,20,50,80,95,100)/100)
x = agg$hours
ys = agg$pause
q99 <- function(x, ...) { x <- quantile(x,probs=c(0.99)) }
v99 = q99(gcdata$Pause.s.)
vmed = median(gcdata$Pause.s.)
plotAreaCorridor(x, ys,ylim=c(0,v99*1.5))
abline(h=vmed, col="lightblue")
abline(h=v99, col="grey")
label=paste("99%=",round(v99,digits=3),"s n=", length(gcdata$date),sep="")
text(x=30, y=v99, label, col="grey", pos=3)
title("NewPar Collection Activity")