2013-08-23 93 views
1

如果數據大於某個值,我想在圖中寫一些條件統計量。ggplot的條件stat_summary R

與傑克瑞安的熱心幫助(Cut data and access groups to draw percentile lines),我可以創建以下腳本組數據到小時,並繪製結果:

# Read example data 
A <- read.csv(url('http://people.ee.ethz.ch/~hoferr/download/data-20130812.csv')) 

# Libraries 
library(doBy) 
library(ggplot2) 
library(plyr) 
library(reshape2) 
library(MASS) 
library(scales) 

# Sample size function 
give.n <- function(x){ 
     return(c(y = min(x) - 0.2, label = length(x))) 
} 

# Calculate gaps 
gaps <- rep(NA, length(A$Timestamp)) 
times <- A$Timestamp 
loss <- A$pingLoss 
gap.start <- 1 
gap.end <- 1 
for(i in 2:length(A$Timestamp)) 
{ #For all rows 
    if(is.na(A$pingRTT.ms.[i])) 
    { #Currently no connection 
     if(!is.na(A$pingRTT.ms.[i-1])) 
     { #Connection lost now 
      gap.start <- i 
     } 
     if(!is.na(A$pingRTT.ms.[i+1])) 
     { # Connection restores next time 
      gap.end <- i+1 
      gaps[gap.start] <- as.numeric(A$Timestamp[gap.end]-A$Timestamp[gap.start], units="secs") 
      loss[gap.start] <- gap.end - gap.start 
     } 
    }  
}    
H <- data.frame(times, gaps, loss) 
H <- H[complete.cases(H),] 
C <- H  
C$dates <- strptime(C$times, "%Y-%m-%d %H:%M:%S") 
C$h1 <- C$dates$hour 

# Calculate percentiles 
cuts <- c(1, .75, .5, .25, 0) 
c <- ddply(C, .(h1), function (x) { summarise(x, y = quantile(x$gaps, cuts)) }) 
c$cuts <- cuts 
c <- dcast(c, h1 ~ cuts, value.var = "y") 
c.melt <- melt(c, id.vars = "h1") 

p <- ggplot(c.h1.melt, aes(x = h1, y = value, color = variable)) + 
geom_point(size = 4) + 
stat_summary(fun.data = max.n, geom = "text", fun.y = max, colour = "red", angle = 90, size=4) + 
scale_colour_brewer(palette="RdYlBu", name="Percentile", guide = guide_legend(reverse=TRUE)) + 
scale_x_continuous(breaks=0:23, limits = c(0,23)) + 
annotation_logticks(sides = "lr") + 
theme_bw() + 
scale_y_log10(breaks=c(1e0,1e1,1e2,1e3,1e4), labels = trans_format("log10", math_format(10^.x)), limits=c(1e0,1e4)) + 
xlab("Hour of day") + ylab("Ping gaps [s]") 
p 

p <- ggplot(c.m1.melt, aes(x = m1/60, y = value, color = variable)) + 
geom_point(size = 1) + 
stat_summary(fun.data = give.n, geom = "text", fun.y = median, angle = 90, size=4) + 
stat_summary(fun.data = max.n, geom = "text", fun.y = max, colour = "red", angle = 90, size=4) + 
scale_colour_brewer(palette="RdYlBu", name="Percentile", guide = guide_legend(reverse=TRUE)) + 
scale_x_continuous(breaks=0:23, limits = c(0,24)) + 
annotation_logticks(sides = "lr") + 
theme_bw() + 
scale_y_log10(breaks=c(1e0,1e1,1e2,1e3,1e4), labels = trans_format("log10", math_format(10^.x)), limits=c(1e0,1e4)) + 
xlab("Time of day") + ylab("Ping gaps [s]") 
p 

這就造成差距與長度每小時分組情節最長的差距書面旁邊的數據點:

Hourly grouped plot without sample number per group

下面是微小的分組情節。這個數字是不可讀的,爲什麼我想添加條件統計數據,如果差距超過5分鐘,或者只有10個最長的差距或類似的東西。

Minutely grouped plot with unreadable stats

我試着只是改變統計功能

max.n.filt <- function(x){ 
    filter = 300 
    if (x > filter) { 
     return(c(y = max(x) + 0.4, label = round(max(10^x),2))) 
    } else { 
     return(c(y=x, label = "")) 
    } 
} 

,並使用這個細微的情節分組。但我得到這個錯誤:

Error in list_to_dataframe(res, attr(.data, "split_labels")) : 
    Results do not have equal lengths 
In addition: There were 50 or more warnings (use warnings() to see the first 50) 
Error in if (nrow(layer_data) == 0) return() : argument is of length zero 
Calls: print ... print.ggplot -> ggplot_gtable -> Map -> mapply -> <Anonymous> 
In addition: Warning message: 
Removed 6 rows containing missing values (geom_point). 

此外,在每小時情節,我想寫旁邊的縫隙的長度每小時的樣本數量。我想我可以在c數據框中添加一個新列,但不幸的是我找不到這樣做的方法。

任何幫助,非常感謝。

回答

2

請參閱?stat_summary。

fun.data : Complete summary function. Should take data frame as input and return data frame as output

你的功能max.n.filt使用試圖評估條件x > filterif()聲明。但是當length(x) > 1時,if()語句僅評估第一個值x的條件。當在數據框上使用時,這將返回從原始輸入x拼湊在一起的列表以及if()語句返回的任何標籤。

​​

嘗試使用ifelse()的函數:

max.n.filt2 <- function(x){ 
    filter = 300     # whatever threshold 
    y = ifelse(x > filter, max(x) + 1, x[,1]) 
    label = ifelse(x > filter, round(max(x),2), NA) 
    return(data.frame(y=y[,1], label=label[,1])) 
} 

> max.n.filt2(data.frame(x=c(10,15,400))) 
    y label 
1 10 NA 
2 15 NA 
3 401 400 

或者,你可能發現它更容易使用geom_text()。我無法重現你的榜樣,但這裏有一個模擬數據集:

set.seed(101) 
sim_data <- expand.grid(m1=1:1440, variable=factor(c(0,0.25,0.5,0.75,1))) 
sim_data$sample_size <- sapply(1:1440, function(.) sample(1:25, 1, replace=T)) 
sim_data$value = t(sapply(1:1440, function(.) quantile(rgamma(sim_data$sample_size, 0.9, 0.5),c(0,0.25,0.5,0.75,1))))[1:(1440*5)] 

只需使用subset論點geom_text()選擇那些你想標記點:

ggplot(sim_data, aes(x = m1/60, y = value, color = variable)) + 
geom_point(size = 4) + geom_text(aes(label=round(value)), subset = .(variable == 1 & value > 25), angle = 90, size = 4, colour = "red", hjust = -0.5) 

如果你有一列樣品的尺寸,那些可併入labelpaste()

ggplot(sim_data, aes(x = m1/60, y = value, color = variable)) + 
geom_point(size = 4) + geom_text(aes(label=paste(round(value),", N=",sample_size)), subset = .(variable == 1 & value > 25), angle = 90, size = 4, colour = "red", hjust = -0.25) 

(或創建一個單獨的列的數據與你想要的任何標籤。)如果你問如何取回樣本大小,您可以修改您的來電ddply()這樣的:

... 
c2 <- ddply(C, .(h1), function (x) { cbind(summarise(x, y = quantile(x$gaps, cuts)), n=nrow(x)) }) 
c2$cuts <- cuts 
c2 <- dcast(c2, h1 + n ~ cuts, value.var = "y") 
c2.h1.melt <- melt(c2, id.vars = c("h1","n")) 
...