How to smooth ecdf plots in r

Question

I have a df with 5 variables,

head(df,15)

               junc  N1.ir  N2.ir    W1.ir    W2.ir    W3.ir
1  pos$chr1:3197398  0.000000  0.000000  0.000000  0.000000  0.000000
2  pos$chr1:3207049  0.000000  0.000000  0.000000  0.000000  0.000000
3  pos$chr1:3411982  0.000000  0.000000  0.000000  0.000000  0.000000
4  pos$chr1:4342162  0.000000  0.000000  0.000000  0.000000  0.000000
5  pos$chr1:4342918  0.000000  0.000000  0.000000  0.000000  0.000000
6  pos$chr1:4767729 -4.369234 -5.123382 -4.738768 -4.643856 -5.034646
7  pos$chr1:4772814 -3.841302 -3.891419 -4.025029 -3.643856 -3.184425
8  pos$chr1:4798063 -5.038919 -4.847997 -5.497187 -4.035624 -7.543032
9  pos$chr1:4798567 -4.735325 -5.096862 -3.882643 -3.227069 -4.983808
10 pos$chr1:4818730 -8.366322 -7.118941 -8.280771 -6.629357 -6.876517
11 pos$chr1:4820396 -5.514573 -6.330917 -5.898853 -4.700440 -5.830075
12 pos$chr1:4822462 -5.580662 -6.914883 -5.562242 -5.380822 -5.703211
13 pos$chr1:4827155 -4.333273 -4.600904 -4.133399 -4.012824 -3.708345
14 pos$chr1:4829569 -4.287866 -3.874469 -3.977280 -4.209453 -4.490326
15 pos$chr1:4857613 -6.902074 -6.074141 -6.116864 -3.989946 -6.474259

Few lines after using melt

> head(ir.m)
              junc variable     value
1 pos$chr1:3197398 N1.ir  0.000000
2 pos$chr1:3207049 N1.ir  0.000000
3 pos$chr1:3411982 N1.ir  0.000000
4 pos$chr1:4342162 N1.ir  0.000000
5 pos$chr1:4342918 N1.ir  0.000000
6 pos$chr1:4767729 N1.ir -4.369234

And summary

> summary(ir)
                 junc           N1.ir          N2.ir           W1.ir       
 neg$chr1:100030088:     1   Min.   :-11.962   Min.   :-12.141   Min.   :-11.817  
 neg$chr1:100039873:     1   1st Qu.: -4.379   1st Qu.: -4.217   1st Qu.: -4.158  
 neg$chr1:10023338 :     1   Median : -2.807   Median : -2.663   Median : -2.585  
 neg$chr1:10024088 :     1   Mean   : -2.556   Mean   : -2.434   Mean   : -2.362  
 neg$chr1:10025009 :     1   3rd Qu.:  0.000   3rd Qu.:  0.000   3rd Qu.:  0.000  
 neg$chr1:10027750 :     1   Max.   : 17.708   Max.   : 16.162   Max.   : 16.210  
 (Other)           :113310                                                        
     W2.ir            W3.ir       
 Min.   :-12.194   Min.   :-11.880  
 1st Qu.: -3.078   1st Qu.: -4.087  
 Median : -1.000   Median : -2.711  
 Mean   : -1.577   Mean   : -2.370  
 3rd Qu.:  0.000   3rd Qu.:  0.000  
 Max.   : 17.562   Max.   : 16.711

I am trying to plot cumulative probability using ggplot and stat_ecdf,

using this code

ggplot(ir.m, aes(x=value)) + stat_ecdf(aes(group=variable,colour = variable))

Plot looks like this,

How do I obtain a smooth curve? Do I need to perform more statistical operation to obtain that?

updated code

ir.d = as.data.frame(ir.m)
denss = split(ir.d, ir.d$variable) %>%
  map_df(function(dw) {
    denss = density(dw$value, from=min(ir.d$value) - 0.05*diff(range(ir.d$value)), 
                   to=max(ir.d$value) + 0.05*diff(range(ir.d$value)))
    data.frame(x=denss$x, y=denss$y, cd=cumsum(denss$y)/sum(denss$y), group=dw$variable[1])
    head(denss)
  })
summary(denss)
> summary(denss)
       x                 y                   cd               group    
 Min.   :-13.689   Min.   :0.0000000   Min.   :0.00000   N1.ir:512  
 1st Qu.: -5.466   1st Qu.:0.0000046   1st Qu.:0.07061   N2.ir:512  
 Median :  2.757   Median :0.0002487   Median :0.99552   W1.ir  :512  
 Mean   :  2.757   Mean   :0.0303942   Mean   :0.65315   W2.ir  :512  
 3rd Qu.: 10.980   3rd Qu.:0.0148074   3rd Qu.:0.99997   W3.ir  :512  
 Max.   : 19.203   Max.   :0.9440592   Max.   :1.00000

plot

ggplot() +
  stat_ecdf(data=ir.d, aes(x, colour=variable), alpha=0.8) +
  geom_line(data=denss, aes(x, cd, colour=group)) +
  theme_classic()

Sorry, I have posted a wrong answer. If it is not the graph of some simulation, you cannot improve this graph. — Viktor, Jan 04 '18 at 17:41
That's the whole point of ECDF... Perhaps you need to rephrase your question and specify what is the underlying question you're after. — Roman Luštrik, Jan 04 '18 at 19:10
@RomanLuštrik I am trying to plot curves between 3 groups (N) and 2(W) here, using ecdf — novicebioinforesearcher, Jan 04 '18 at 19:53

eipi10 · Accepted Answer · 2021-11-12T00:38:15.473

13

The ecdf follows the data exactly, without any smoothing. However, you can create a smoothed cumulative density by generating a kernel density estimate (basically a smoothed histogram) from the data and creating an "ecdf" from that. Here's an example with fake data:

First we generate a kernel density estimate using the density function. This gives us, by default, a density estimate on a grid of 512 x-values. Then we use that as the "data" for calculating the ecdf, which is just the cumulative sum of the density (or, for any given point a along the x axis, the value of the ecdf at a is the area under the kernel density curve (that is, the integral from -Inf to a).

I've packaged the code into a function below so you can see how changing the adjust parameter of the density function changes the smoothed ecdf. A smaller value of adjust reduces the amount of smoothing, creating a density estimate that more closely follows the data. You can see in the plots below that setting adj=0.1 results in less smoothing of the smoothed ecdf so that it more closely follows the step in the original ecdf.

library(ggplot2)

smooth_ecd = function(adj = 1) {

  # Fake data
  set.seed(2)       
  dat = data.frame(x=rnorm(15))
  
  # Extend range of density estimate beyond data
  e = 0.3 * diff(range(dat$x))
  
  # Kernel density estimate of fake data
  dens = density(dat$x, adjust=adj, from=min(dat$x)-e, to=max(dat$x) +e)
  dens = data.frame(x=dens$x, y=dens$y)
  
  # Plot kernel density (blue), ecdf (red) and smoothed ecdf (black)
  ggplot(dat, aes(x)) + 
    geom_density(adjust=adj, colour="blue", alpha=0.7) +
    geom_line(data=dens, aes(x=x, y=cumsum(y)/sum(y)), size=0.7, colour='grey30') +
    stat_ecdf(colour="red", size=0.6, alpha=0.6) +
    theme_classic() +
    labs(title=paste0("adj=",adj))
}

smooth_ecd(adj=1)
smooth_ecd(adj=0.3)
smooth_ecd(adj=0.1)

Here's some code for doing this by group:

library(tidyverse)

# Fake data with two groups
set.seed(2)
dat = data.frame(x=c(rnorm(15, 0, 1), rnorm(20, 0.2, 0.8)), 
                 group=rep(LETTERS[1:2], c(15,20)))

# Split the data by group and calculate the smoothed cumulative density for each group
dens = split(dat, dat$group) %>% 
  map_df(function(d) {
    dens = density(d$x, adjust=0.1, from=min(dat$x) - 0.05*diff(range(dat$x)), 
                   to=max(dat$x) + 0.05*diff(range(dat$x)))
    data.frame(x=dens$x, y=dens$y, cd=cumsum(dens$y)/sum(dens$y), group=d$group[1])
  })

Now we can plot each smoothed cumulative density. In the plot below, I've included a call to stat_ecdf with the original data for comparison.

ggplot() +
  stat_ecdf(data=dat, aes(x, colour=group), alpha=0.8, lty="11") +
  geom_line(data=dens, aes(x, cd, colour=group)) +
  theme_classic()

UPDATE: Using your data sample, here's what I get. I have no idea how you got that long nucleotide string as the x-value in your plot, as a variable like that doesn't appear anywhere in the data you posted.

# Melt data
dat = gather(df, variable, value, -junc)

# Split the data by group and calculate the smoothed cumulative density for each group
dens = split(dat, dat$variable) %>% 
  map_df(function(d) {
    dens = density(d$value, adjust=0.1, from=min(dat$value) - 0.05*diff(range(dat$value)), 
                   to=max(dat$value) + 0.05*diff(range(dat$value)))
    data.frame(x=dens$x, y=dens$y, cd=cumsum(dens$y)/sum(dens$y), group=d$variable[1])
  })

ggplot() +
  stat_ecdf(data=dat, aes(value, colour=variable), alpha=0.8, lty="11") +
  geom_line(data=dens, aes(x, cd, colour=group)) +
  theme_classic()

edited Nov 12 '21 at 00:38

answered Jan 04 '18 at 19:31

eipi10

91,525
24
209
285

I would refrain from calling that smooth_ecdf since it is not a function. The ecdf-function returns a function hence its name. Call yours `smooth_plot` or some such. – IRTFM Jan 04 '18 at 19:36
Does `smooth_ecd` give you a sense of inner calm? – eipi10 Jan 04 '18 at 19:39
My (tongue-in-cheek) comment was for @42-. Let me know what you're confused about and I'll see if I can help. – eipi10 Jan 04 '18 at 20:23
error : Error in function_list[[k]](value) : could not find function "map_df" – novicebioinforesearcher Jan 04 '18 at 21:20
1

Make sure you load the tidyverse package (`library(tidyverse)`). – eipi10 Jan 04 '18 at 21:23
there were three plots right did you delete something? – novicebioinforesearcher Jan 04 '18 at 21:32
1

There are still three plots. Two together in the middle and one at the end. – eipi10 Jan 04 '18 at 21:32
No the one with smooth lines – novicebioinforesearcher Jan 04 '18 at 21:33
I haven't deleted anything from my answer. All of the plots have lines for the smoothed versions of the ecdf. You can look at the [edit chain](https://stackoverflow.com/posts/48102268/revisions) to see all revisions to the answer. – eipi10 Jan 04 '18 at 21:34
Ahh thanks!apologies for the confusion, `0.05` a constant or was it based on something? – novicebioinforesearcher Jan 04 '18 at 21:37
When i try it on my data i get this error `Error in data.frame(x = dens$value, y = dens$y, cd = cumsum(dens$y)/sum(dens$y), : arguments imply differing number of rows: 0, 512, 1 ` – novicebioinforesearcher Jan 04 '18 at 21:44
1

`0.05` is just to extend the x-range by 5% in each direction so that the smoothed ecdf plot includes the full range of the data plus a bit extra. – eipi10 Jan 04 '18 at 21:49
1

As for the error, it looks like `dens$value` is the problem. It should be `dens$x`. Note that this section of the code is using the `x` and `y` values returned by the `density` function, not a column from your data frame. – eipi10 Jan 04 '18 at 21:51
That did it, but the plot has something weird on x axis, please see edit – novicebioinforesearcher Jan 04 '18 at 22:11
It looks like the x-axis in your second plot is a categorical variable (it looks like a very long nucleotide chain). Hard to be sure what's going on without a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). – eipi10 Jan 04 '18 at 22:23
no way i can attach the original file here(file size), is there a way i can send it across to you. – novicebioinforesearcher Jan 04 '18 at 22:27
1

We just need a few rows for each of two different groups. Use `dput`. For example, paste into your question the output of `dput(ir.m[c(1:5, 21:25), ])` to provide rows 1 through 5 and rows 21 through 25. Just pick a few rows that represent two different levels of `variable`. – eipi10 Jan 04 '18 at 22:30
@eipi10 done, i pasted the original data(without melt) 15 rows. – novicebioinforesearcher Jan 04 '18 at 22:39
accepted it, i have no idea as well why i get that string, i did dev.off() now its gone – novicebioinforesearcher Jan 04 '18 at 23:30

score -2 · Answer 2 · answered Jun 18 '20 at 14:26

-2

This is an older thread, however, I just want to mention that stat_ecdf(..., geom = "line") might be a suitable solution for some people to avoid the steps from geom_step in the ecdf curve. -Michael

answered Jun 18 '20 at 14:26

mod_che

137
1
6

How to smooth ecdf plots in r

updated code

plot

2 Answers2