Calculate probability_trans(distribution="norm") without ggplot

Question

I am having trouble understanding an R plotting script one of my coworkers provided me. I understand that ggplot can be used to transform data in the plotting routine, but I am looking for an alternate method to calculate the transform using other functions in r. For the data bwpeaks

structure(list(PROB = c(0.25, 0.116666666666667, 0.15, 0.366666666666667, 
0.4, 0.133333333333333, 0.633333333333333, 0.5, 0.483333333333333, 
0.516666666666667, 0.1, 0.3, 0.666666666666667, 0.9, 0.716666666666667, 
0.466666666666667, 0.283333333333333, 0.583333333333333, 0.0833333333333333, 
0.35, 0.416666666666667, 0.866666666666667, 0.6, 0.2, 0.75, 0.533333333333333, 
0.933333333333333, 0.566666666666667, 0.816666666666667, 0.966666666666667, 
0.433333333333333, 0.85, 0.95, 0.333333333333333, 0.8, 0.0666666666666667, 
0.733333333333333, 0.883333333333333, 0.183333333333333, 0.0166666666666667, 
0.7, 0.45, 0.266666666666667, 0.216666666666667, 0.616666666666667, 
0.166666666666667, 0.916666666666667, 0.383333333333333, 0.683333333333333, 
0.316666666666667, 0.05, 0.766666666666667, 0.233333333333333, 
0.0333333333333333, 0.983333333333333, 0.783333333333333, 0.833333333333333, 
0.65, 0.55), FLOW = c(2570, 1330, 1500, 3810, 4100, 1400, 10900, 
5590, 5520, 6460, 1240, 2960, 11100, 33400, 13100, 4710, 2730, 
9640, 1200, 3330, 4130, 24500, 9650, 2260, 13900, 6980, 38800, 
8300, 17700, 49600, 4170, 22200, 47700, 3310, 15900, 1020, 13500, 
25200, 2240, 135, 12500, 4650, 2580, 2300, 9680, 1660, 34600, 
4010, 11800, 3070, 675, 14600, 2400, 406, 64300, 15100, 19800, 
10900, 8040)), class = "data.frame", row.names = c(NA, -59L))

where bwpeaks$PROB is the weibull plotting position for the vector bwpeaks$FLOW. The ploting routine looks like:


library(ggplot2)
library(scales)
library(lmomco)

log10_minor_break = function (...){
  function(x) {
    minx         = floor(min(log10(x), na.rm=T))-1;
    maxx         = ceiling(max(log10(x), na.rm=T))+1;
    n_major      = maxx-minx+1;
    major_breaks = seq(minx, maxx, by=1)
    minor_breaks = 
      rep(log10(seq(1, 9, by=1)), times = n_major)+
      rep(major_breaks, each = 9)
    return(10^(minor_breaks))
  }
}
  xbreaks <- c(0.002,0.01,0.10,0.25,0.5,0.8,0.9,0.95,0.975,0.99,0.995, 0.998)
  log.range <- log10(range(bwpeaks$FLOW, na.rm = TRUE))
  lower <- 10^floor(log.range[1])
  upper <- 10^ceiling(log.range[2])
  cap <- 100
  ybreaks <- NULL
  while(cap < upper) {
    ybreaks <- c(ybreaks, seq(cap, cap*1, by = cap))
    cap <- cap * 10
  }
  ggplot(bwpeaks) + 
    geom_point(aes(x=PROB, y=FLOW)) + 
    theme_bw() + 
    scale_y_continuous(trans="log10", 
                       breaks=round(ybreaks, digits=2), 
                       minor_breaks=log10_minor_break(), 
                       name="Discharge (cfs)", labels=comma) +
    scale_x_continuous(trans=probability_trans(distribution="norm"),
                       breaks=xbreaks, 
                       labels=signif(prob2T(xbreaks), digits=3),
                       name="Return period (yrs)")

I am confused by scale_x_continuous(trans=probability_trans(distribution="norm"). Is there a way to calculate the transformed x variable bwpeaks$PROB without using ggplot?

user_8583393 · Answer

The function probability_trans comes from the scales package. It returns a list.

tmp <- scales::probability_trans(distribution="norm")
str(tmp)
#List of 7
# $ name        : chr "prob-norm"
# $ transform   :function (x)  
# $ inverse     :function (x)  
# $ breaks      :function (x, n = n_default)  
# $ minor_breaks:function (b, limits, n)  
# $ format      :function (x)  
# $ domain      : num [1:2] -Inf Inf
# - attr(*, "class")= chr "trans"

So you could try

tmp$transform(bwpeaks$PROB)
#[1] -0.67448975 -1.19181617 -1.03643339 -0.34069483 -0.25334710
#[6] -1.11077162  0.34069483  0.00000000 -0.04178930  0.04178930
# ...

user_13095326 · Answer

@markus's explanation is great. If you read the source code you'll see that it calls match.fun for paste0("q",distribution). So you're looking for the qnorm function.

probability_trans(distribution="norm")$transform(bwpeaks$PROB)
 [1] -0.67448975 -1.19181617 -1.03643339 -0.34069483 -0.25334710 -1.11077162  0.34069483  0.00000000 -0.04178930  0.04178930 -1.28155157 -0.52440051
[13]  0.43072730  1.28155157  0.57296755 -0.08365173 -0.57296755  0.21042839 -1.38299413 -0.38532047 -0.21042839  1.11077162  0.25334710 -0.84162123
[25]  0.67448975  0.08365173  1.50108595  0.16789400  0.90273479  1.83391464 -0.16789400  1.03643339  1.64485363 -0.43072730  0.84162123 -1.50108595
[37]  0.62292572  1.19181617 -0.90273479 -2.12804523  0.52440051 -0.12566135 -0.62292572 -0.78350038  0.29673784 -0.96742157  1.38299413 -0.29673784
[49]  0.47704043 -0.47704043 -1.64485363  0.72791329 -0.72791329 -1.83391464  2.12804523  0.78350038  0.96742157  0.38532047  0.12566135

qnorm(bwpeaks$PROB)
 [1] -0.67448975 -1.19181617 -1.03643339 -0.34069483 -0.25334710 -1.11077162  0.34069483  0.00000000 -0.04178930  0.04178930 -1.28155157 -0.52440051
[13]  0.43072730  1.28155157  0.57296755 -0.08365173 -0.57296755  0.21042839 -1.38299413 -0.38532047 -0.21042839  1.11077162  0.25334710 -0.84162123
[25]  0.67448975  0.08365173  1.50108595  0.16789400  0.90273479  1.83391464 -0.16789400  1.03643339  1.64485363 -0.43072730  0.84162123 -1.50108595
[37]  0.62292572  1.19181617 -0.90273479 -2.12804523  0.52440051 -0.12566135 -0.62292572 -0.78350038  0.29673784 -0.96742157  1.38299413 -0.29673784
[49]  0.47704043 -0.47704043 -1.64485363  0.72791329 -0.72791329 -1.83391464  2.12804523  0.78350038  0.96742157  0.38532047  0.12566135

Calculate probability_trans(distribution="norm") without ggplot

2 Answers