3
$\begingroup$

I have df1

> df1
   V1      V2      V3      V4   V5     V6       V7  V8
1   1   EFNB2   EFNB2   EFNB2   72  0.667  0.05980 YES
2   2    AOC2    AOC2    AOC2  108  0.609  0.13000 YES
3   3     FOS     FOS     FOS  142  0.554  0.19300 YES
4   4   CCNE1   CCNE1   CCNE1  167  0.530  0.25700 YES
5   5    CBX4    CBX4    CBX4  275  0.448  0.26700 YES
6   6    CRY1    CRY1    CRY1  343  0.400  0.29000 YES
7   7 TNFAIP6 TNFAIP6 TNFAIP6  600  0.194  0.18500  NO
8   8    SSX3    SSX3    SSX3  641  0.171  0.18900  NO
9   9    WNT2    WNT2    WNT2 1211 -0.173 -0.08030  NO
10 10   FGF18   FGF18   FGF18 1296 -0.225 -0.09130  NO
11 11    PIM2    PIM2    PIM2 1386 -0.273 -0.09780  NO
12 12   PLCL1   PLCL1   PLCL1 1636 -0.415 -0.16700  NO
13 13     MX2     MX2     MX2 1676 -0.437 -0.12400  NO
14 14   ABCA6   ABCA6   ABCA6 1720 -0.457 -0.07980  NO
15 15 HSPA12A HSPA12A HSPA12A 1788 -0.501 -0.04190  NO
16 16 TNFSF10 TNFSF10 TNFSF10 1944 -0.855  0.00155  NO

I am supposed to iterate through RES column and find the point of inflection and mark everything from there with YES -CORE_ENRICHMENT-. P-value is identified and according to my very little knowledge in statistics that we need to set the second derivative of f(x) to 0 in order to identify it. Is there anything that I am missing or a formula that might help me to understand how can I identify it? I need to implement a program in R while I am really unable to understand the problem.

Edit: Based on User603 suggestion: I modified the code as follow, and plotted this plot:

df1 <- read.table("file.XLS", header = FALSE)
mlearn <- df1$V7
#plot of RES:
plot(mlearn, ylim = c(-0.2, 0.3), type = 'p', pch = 16)
#plot of RES^1:
#plot of the smoothed values of RES:
lines(fdata.deriv(mlearn, nderiv = 0, method = "bspline", class.out = 'fd', nbasis = 6), col = 'blue')
#the first derivative of the smoothed values of RES
lines(fdata.deriv(mlearn, nderiv = 1, method = "bspline", class.out = 'fd', nbasis = 6), col = 'red') 
legend('topleft', col = c('blue', 'red'), lty = c(1, 1), legend = c('Smoothed curve', 'First derivative of smoothed curve'))
abline(h = 0)

enter image description here

1- If I were to identify IP from the graph, it should be something around 5, however, I've been told by my doctor that the other crossing point around 12 should be an IP as well.

I used this code to check the condition: All redline points cross 0.

unlins <- RedLine_points$data[1, 1:16]
ifelse(diff(sign(unlins))==-2|diff(sign(unlins))==2 , mlearn, NA)

However, when I run the code, what returns is two points: First point is the point before the IP around 5 which I believe is wrong. 2nd point is exactly the point around 12.

How can I enhance my code to return the exact IPs. I need a general formula as I am working on 18 other files.

Update: This is the modified code suggested

install.packages("matrixStats")
install.packages('fda.usc')
library(matrixStats)
library(fda.usc)
df1 <- read.table("file.XLS", header = FALSE)
mlearn <- df1$V7
rc_container <- fdata.deriv(mlearn, nderiv = 1, method = "bspline", class.out = 'fdata', nbasis = 6)
red_curve    <- cbind(rc_container$argvals, t(rc_container$data))
red_curve    <- data.frame(red_curve)
colnames(red_curve) <- c('x', 'y');
plot(red_curve$x, red_curve$y * 30, type = 'l', col = 'red')
diff_rc <- cbind(2:length(red_curve$y)-1 , 
                 red_curve$y[-c(length(red_curve$y))],
                 red_curve$y[-c(1)]) 
x_loc <- red_curve$x[diff_rc[which(rowProds(diff_rc[, 2:3]) < 0), 1]] #return index "index"
y_loc <- red_curve$y[diff_rc[which(rowProds(diff_rc[, 2:3]) < 0), 1]]
abline(v = x_loc)
abline(h = 0)

And this is the result

> diff_rc
      [,1]         [,2]         [,3]
 [1,]    1  0.043131541  0.079957068
 [2,]    2  0.079957068  0.074712332
 [3,]    3  0.074712332  0.041978947
 [4,]    4  0.041978947 -0.003661470
 [5,]    5 -0.003661470 -0.047627302
 [6,]    6 -0.047627302 -0.077841351
 [7,]    7 -0.077841351 -0.092244072
 [8,]    8 -0.092244072 -0.091280341
 [9,]    9 -0.091280341 -0.075395031
[10,]   10 -0.075395031 -0.045033013
[11,]   11 -0.045033013 -0.003151946
[12,]   12 -0.003151946  0.037239385
[13,]   13  0.037239385  0.060619410
[14,]   14  0.060619410  0.051466559
[15,]   15  0.051466559 -0.005740738
> x_loc
[1]  4 12 15
> y_loc
[1]  0.041978947 -0.003151946  0.051466559

enter image description here

Can I also use colProds to check sign changes of diff_rc[,2] and return the minimum value? My issue with rawProds is that it doesn't necessarily return the same point. It returns where the change of the sign has occurred but not where the point is exactly at.

Example: Based on the example above, x_loc 4 is where the change of sign occurred, while what I need is x_loc = 5 because y_loc= -0.003661470 is the fifth point and it is the one crosses 0 not 0.041978947.

When I try to use colProds, I get an error:

> colProds(diff_rc[,2])
Error in double(length = n) : invalid 'length' argument
In addition: Warning message:
Argument 'x' is of class ‘numeric’, but should be a matrix. The use of a ‘numeric’ is not supported, the correctness of the result is not guaranteed, and will be defunct (produce an error) in a future version of matrixStats. Please update your code according

Output of this method and the one suggested here

> Identify_IP(df)
[[1]]
[1]  0.2570 -0.1670 -0.0419

[[2]]
[1]  4 12 15

This is the new indexes

> x_loc
[1]  4.5 12.5 15.5
> inflection_points <- mlearn[x_loc]
> inflection_points
[1]  0.2570 -0.1670 -0.0419

X_loc has changed but inflection points of both methods are the same.

$\endgroup$
12
  • $\begingroup$ I am back on this problem as I was having my exams. Can you tell me if the condition I added is correct? Or there is another condition that can be used? 2- Can the point around 12 be an IP as well according to my doctor's? $\endgroup$
    – MarJamil
    Oct 29, 2017 at 6:56
  • $\begingroup$ Your doctor is right. ~12 is also an inflection point. But ~12 is not a local maximum, it is a local minimum. If you want all inflection points, then you should select all places where the sign of the red curve changes (either from positive to negative or vis-versa). $\endgroup$
    – user603
    Oct 30, 2017 at 14:33
  • $\begingroup$ I really appreciate your help. I apologize for the too much questions. You are doing a great help to someone who knows nothing about R and very little about statistics. I really made a progress out of your answers. Can't thank you enough :) $\endgroup$
    – MarJamil
    Nov 2, 2017 at 12:12
  • $\begingroup$ Any updates on this? I have been trying different approaches but nothing has worked so far. $\endgroup$
    – MarJamil
    Nov 6, 2017 at 9:16
  • $\begingroup$ OK, I will rewrite the code for your data. Which one of the column in your data is RES? $\endgroup$
    – user603
    Nov 6, 2017 at 9:22

1 Answer 1

11
$\begingroup$

The problem with the approach the OP proposes is that it fails to acknowledge the random nature of the data. In other words, that the value of RES are only an observed sample, not the whole data.

The inflection point computed in this way and based on a sample of values of $\text{RES}$ -- call it $\text{IP}(\text{RES})$-- is merely an estimate of the true inflection point --call it $\text{IP}_0$--: the one we would have obtained had we had access to all the data (not just one sample).

See, if we had drawn another sample $\text{RES}^1$ from the same distribution and used that one as data instead, we would obtain a different value of $\text{IP}(\text{RES}^1)$ too.

The insight is to recognize that $\text{IP}(\text{RES})$ is not necessarily the best estimate of $\text{IP}_0$. It is typically too noisy for that. Noisy here means that it is geared towards describing artifacts of $\text{RES}$ that may not be present in another sample $\text{RES}^1$.

Typically, one smooths the value of the variable of interest (here $\text{RES}$) with a scatter-plot smoother (there are many out there but we must choose one with explicit derivatives like bsplines or Fourrier).

Then, you just look for locations where the first derivative of the smoothed curve is 0 (and for a local maximum, the first derivative of the smoothed curve is positive to the left and negative to the right of that location). This yields smoother -- less noisy-- estimate of the true inflection point $\text{IP}_0$.

Here is some R illustrated example (to keep the code readable, I chose to defer the actual computation of the derivative of the scatterplot smoother to the fda.usc::fdata.deriv function but it could also be done by hand);

library(fda.usc)
data(phoneme)
n = length(phoneme$learn[1, ]$data)

set.seed(123)
test_sample <- sample.int(n, n/2)
fit_sample  <- (1:n)[-test_sample]
mlearn      <- phoneme$learn[1, fit_sample]  #think of this as RES
mtest       <- phoneme$learn[1, test_sample] #think of this as RES^1
#plot of RES:
plot(mlearn, ylim = c(-5, 20), type = 'p', pch = 16)
#plot of RES^1:
points(mtest$argvals, mtest$data, pch = 8)
#plot of the smoothed values of RES:
lines(fdata.deriv(mlearn, nderiv = 0, method = "bspline", class.out = 'fd', nbasis = 6), col = 'blue')
#plot of 30 times (so it is easily visible) 
#the first derivative of the smoothed values of RES
lines(fdata.deriv(mlearn, nderiv = 1, method = "bspline", class.out = 'fd', nbasis = 6) * 30, col = 'red') 
#plot of the smoothed values of RES^1:
lines(fdata.deriv(mtest, nderiv = 0, method = "bspline", class.out = 'fd', nbasis = 6), col = 'blue')
#the first derivative of the smoothed values of RES^1
lines(fdata.deriv(mtest, nderiv = 1, method = "bspline", class.out = 'fd', nbasis = 6) * 30, col = 'red') 
legend('topleft', col = c('blue', 'red'), lty = c(1, 1), legend = c('Smoothed curve', 'First derivative of smoothed curve'))
abline(h = 0)

enter image description here

So in this example the original data ($\text{RES}$) are the dots (ordered by frequencies). The stars show another sample of the same distribution ($\text{RES}^1$). The blue line is the result of applying a scatter plot smoother to the entries of ($\text{RES}$) (the dotted blue line is the same but based on $\text{RES}^1$).

The red curve is the (explicit) first derivative of the blue curve (multiplied by 30 so it is visible on the same scale). Here too, the dotted red line is constructed like the red one but using $\text{RES}^1$ instead of $\text{RES}$.

On the red line, there is a single localization where the red curve crosses 0 and is positive to the left and negative to the right of that point (around frequency = 100). This will be the estimate of $\text{IP}_0$ based on $\text{RES}$ ($\text{IP}(\text{RES})$).

You can see that this estimate is close to the maximum values of $\text{RES}$ but doesn't quiet match it. This is because the blue line tries to dampen (or smooth) the idiosyncratic features of ($\text{RES}$).

Because of this dampening, the estimate of $\text{IP}_0$ you would have obtained from $\text{RES}^1$ is not that different from the one you obtained using $\text{RES}$. In this case $\text{IP}(\text{RES}^1)$ is just a little to the right of $\text{IP}(\text{RES})$.

Your estimate of the true inflection point you get in this way is less moved around by the idiosyncrasies of the individual samples you use to compute it and will be a more reliable predictor of the true $\text{IP}_0$.

Edit:

Code to find the points where the red curve crosses 0. This is an adaptation of the code posted by @josliber here.

library(matrixStats)
rc_container <- fdata.deriv(mlearn, nderiv = 1, method = "bspline", class.out = 'fdata', nbasis = 6)
red_curve    <- cbind(rc_container$argvals, t(rc_container$data))
red_curve    <- data.frame(red_curve)
colnames(red_curve) <- c('x', 'y');
plot(red_curve$x, red_curve$y * 30, type = 'l', col = 'red')
 pos <- head(cumsum(rle(red_curve$y >= 0)$lengths), -1)
x_loc <- (red_curve$x[pos] + red_curve$x[pos+1]) / 2
c("positive", "negative")[head(rle(red_curve$y >= 0)$values + 1, -1)]
abline(v = x_loc)
abline(h = 0)

enter image description here

The script above finds all interval where the first derivative of the blue curve (the red curve) crosses 0 (since you now want both local minimae and local maximae).

Edit2:

I am now adapting the code to run on the OP's specific example. According to the OP, the variable of interest is 'V7'.

library(fda.usc)
a1     <- read.table('op_data.txt', header = FALSE)
mlearn <- a1[, 'V7']
plot(mlearn, type = 'p', pch = 16)
lines(fdata.deriv(mlearn, nderiv = 0, method = "bspline", class.out = 'fd', nbasis = 6), col = 'blue')
lines(fdata.deriv(mlearn, nderiv = 1, method = "bspline", class.out = 'fd', nbasis = 6), col = 'red', ylab = 'First derivative of smoothed mlearn') 
abline(h = 0)

enter image description here

Now, it remains to find the points where the red curve crosses 0. Using the code above, we get:

rc_container <- fdata.deriv(mlearn, nderiv = 1, method = "bspline", class.out = 'fdata', nbasis = 6)
red_curve    <- cbind(rc_container$argvals, t(rc_container$data))
red_curve    <- data.frame(red_curve)
colnames(red_curve) <- c('x', 'y');
plot(red_curve$x, red_curve$y, type = 'l', col = 'red')
pos <- head(cumsum(rle(red_curve$y >= 0)$lengths), -1)
x_loc <- (red_curve$x[pos] + red_curve$x[pos+1]) / 2
c("negative", "positive")[head(rle(red_curve$y >= 0)$values + 1, -1)]

enter image description here

So the localizations of changes in sign in the first derivative of the smoothed curve are:

x_loc

which in your case are $x=\{4.5, 12.5, 15.5\}$. If we look at the signs of the red curve between these intervals:

c("negative", "positive")[head(rle(red_curve$y >= 0)$values + 1, -1)]

we see that the blue line is increasing between indexes 1 and 4.5 then decreasing between indexes 4.5 and 12.5 and then increasing again between indexes 12.5 and 16.

$\endgroup$
6
  • 1
    $\begingroup$ Thank you so much for the fascinating explanation. I wasn't able to find a similar answer on the internet for days. I think that I pretty much understand your explanation. I am trying to study the R code you provided and see how can I implement a similar one to my problem. I might need your help here again. Thank you! $\endgroup$
    – MarJamil
    Oct 13, 2017 at 21:05
  • $\begingroup$ No problem, all the functions I used are open source and documented and easy to use and install, see the fda.usc package. There are some ancillary aspects I left out to make the answer more readable, but you should be able to find them addressed on cross validated ('how to select the parameters of a bspline?'). All in all this is fairly standard stuff;) $\endgroup$
    – user603
    Oct 14, 2017 at 14:52
  • $\begingroup$ I have modified my question above @user603 I only used RES samples for now as I don't have the whole one. I'm still a bit confused $\endgroup$
    – MarJamil
    Oct 15, 2017 at 8:35
  • $\begingroup$ I appreciate your help. now, I should check the change of signs for all y_loc (prev and next)points. Since there was 5 points returned of y_loc and I need to make sure which makes IP. Is that correct? $\endgroup$
    – MarJamil
    Oct 30, 2017 at 12:10
  • $\begingroup$ @MarJamil: I updated my code. We only need to consider pairs of succesive points where the sign of the first is not the same as the sign of the second. The changes are in the '#Edit' part only. Check if it helps. Can you update your question to use the new code? $\endgroup$
    – user603
    Oct 30, 2017 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.