Q:

Finding peak value in a bell shaped curve signal using R

I have data such as this. This is for one participant steering for 2.5s. I'm interested in the yaw rate (diff_YR) signal
df <- structure(list(ppid_trialn = c("1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39", 
"1_5_39", "1_5_39", "1_5_39", "1_5_39", "1_5_39"), diff_YR = c(0, 
0, 0, 0, 0, 0, 0, 0, 0, -0.0029370662937063, -0.00220279636363637, 
-0.0015501090909091, 0.00775058834498834, 0.00864802237762236, 
0.0260256428904429, 0.0375291375291375, 0.0518065249417249, 0.0751398564102564, 
0.0939044233100233, 0.108508151048951, 0.134778553939394, 0.149545457902098, 
0.164393939393939, 0.169452214452214, 0.180058275058275, 0.166678321678322, 
0.167494172494173, 0.169207459207459, 0.162272727272727, 0.157051282051282, 
0.158846153846154, 0.157377622377622, 0.14979020979021, 0.14513986013986, 
0.131025641025641, 0.112016317016317, 0.1104662004662, 0.0927622377622378, 
0.0775874125874126, 0.0607808857808858, 0.0456060606060607, 0.0353263403263403, 
0.0183566433566434, 0.00962703962703965, 0.00122377622377622, 
-0.00228438228438226, 0.00367132867132869, -0.00367132867132869, 
-0.0035897435897436, 0.000244755244755215, -0.00832167832167834, 
-0.0115034965034965, -0.012972027972028, -0.0214568764568764, 
-0.0257808857808857, -0.0317365967365967, -0.0406293706293706, 
-0.0457692307692308, -0.0470745920745922, -0.0545804195804196, 
-0.0703263403263403, -0.0793822843822844, -0.0917832167832168, 
-0.109568764568765, -0.117808857808858, -0.131025641025641, -0.143834498834499, 
-0.146363636363636, -0.150198135198135, -0.164801873193473, -0.171573416083916, 
-0.179242416083916, -0.186993001165501, -0.187400928904429, -0.177529136363636, 
-0.171410259254079, -0.171981353379953, -0.160151519347319, -0.153706300699301, 
-0.147587422377622, -0.152156168531469, -0.153379953613054, -0.158682984382284, 
-0.16357808974359, -0.157703964335664, -0.161701633799534, -0.174755243076923, 
-0.171655011655012, -0.178100233100233, -0.17997668997669, -0.190501165501166, 
-0.195559440559441, -0.206002331002331, -0.223135198135198, -0.215874125874126, 
-0.229825174825175, -0.23472027972028, -0.230559440559441, -0.232762237762238, 
-0.225174825174825, -0.212365967365967, -0.196783216783217, -0.18462703962704, 
-0.168881118881119, -0.132983682983683, -0.109160839160839, -0.0757925407925408, 
-0.0495221445221446, -0.0265151515151514, 0.00155011655011648, 
0.0209673659673659, 0.0438111888111889, 0.0527039627039627, 0.0650233100233101, 
0.079055944055944, 0.0835431235431236, 0.0921911421911422, 0.10467365967366, 
0.110547785547786, 0.119277389277389, 0.128414918414918, 0.14, 
0.148648018648019, 0.161456876456876, 0.167167832167832, 0.182424242424242, 
0.201841491841492), timestamp_zero = c(0.300099999999986, 0.316699999999997, 
0.33359999999999, 0.350200000000001, 0.36669999999998, 0.383299999999991, 
0.400000000000006, 0.416599999999988, 0.433300000000003, 0.449999999999989, 
0.466800000000006, 0.483299999999986, 0.5, 0.516699999999986, 
0.5334, 0.550099999999986, 0.566699999999997, 0.583399999999983, 
0.599999999999994, 0.616600000000005, 0.633299999999991, 0.650099999999981, 
0.666599999999988, 0.683300000000003, 0.699999999999989, 0.7166, 
0.733299999999986, 0.75, 0.766699999999986, 0.7834, 0.800099999999986, 
0.816699999999997, 0.833499999999987, 0.849999999999994, 0.86669999999998, 
0.883399999999995, 0.900000000000006, 0.916699999999992, 0.933300000000003, 
0.949999999999989, 0.966700000000003, 0.983299999999986, 1, 1.01669999999999, 
1.0334, 1.04999999999998, 1.0667, 1.08329999999998, 1.1001, 1.11660000000001, 
1.13329999999999, 1.15000000000001, 1.16659999999999, 1.1833, 
1.19999999999999, 1.2166, 1.23319999999998, 1.25, 1.26669999999999, 
1.2833, 1.29999999999998, 1.31659999999999, 1.33329999999998, 
1.34999999999999, 1.36669999999998, 1.38339999999999, 1.40000000000001, 
1.41669999999999, 1.4333, 1.45009999999999, 1.4667, 1.48329999999999, 
1.5, 1.51669999999999, 1.5334, 1.54999999999998, 1.56659999999999, 
1.58329999999998, 1.59999999999999, 1.61660000000001, 1.63339999999999, 
1.65000000000001, 1.66659999999999, 1.6833, 1.69999999999999, 
1.7166, 1.73329999999999, 1.7501, 1.76659999999998, 1.7834, 1.79999999999998, 
1.81659999999999, 1.83329999999998, 1.8501, 1.86660000000001, 
1.88319999999999, 1.8999, 1.91659999999999, 1.9333, 1.94989999999999, 
1.9667, 1.98339999999999, 1.9999, 2.01659999999998, 2.03319999999999, 
2.04990000000001, 2.06659999999999, 2.08320000000001, 2.09989999999999, 
2.11660000000001, 2.13319999999999, 2.15000000000001, 2.16649999999998, 
2.1832, 2.19979999999998, 2.2165, 2.23310000000001, 2.2499, 2.26659999999998, 
2.28319999999999, 2.29999999999998, 2.31649999999999, 2.33320000000001, 
2.34989999999999, 2.3665, 2.38319999999999, 2.3998), SWA = c(0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0900001, 0.1800001, 
0.36, 0.45, 0.81, 0.9900001, 1.35, 1.8, 2.16, 2.61, 3.15, 3.51, 
3.96, 4.41, 4.77, 5.13, 5.67, 6.03, 6.39, 6.75, 7.11, 7.47, 7.65, 
7.92, 8.1, 8.28, 8.46, 8.46, 8.55, 8.55, 8.55, 8.55, 8.55, 8.55, 
8.55, 8.55, 8.55, 8.46, 8.46, 8.46, 8.28, 8.28, 8.19, 8.01, 7.92, 
7.74, 7.56, 7.47, 7.2, 6.84, 6.57, 6.21, 5.85, 5.49, 5.13, 4.68, 
4.23, 3.78, 3.33, 2.79, 2.34, 1.89, 1.44, 1.08, 0.72, 0.27, -0.1799999, 
-0.45, -0.9, -1.26, -1.8, -2.16, -2.61, -3.06, -3.42, -3.96, 
-4.5, -4.95, -5.49, -6.03, -6.57, -7.2, -7.83, -8.46, -8.91, 
-9.63, -10.08, -10.62, -11.16, -11.52, -11.88, -12.15, -12.42, 
-12.51, -12.51, -12.51, -12.51, -12.42, -12.24, -12.06, -11.88, 
-11.61, -11.43, -11.16, -10.89, -10.62, -10.17, -9.81, -9.54, 
-9.09, -8.64, -8.19, -7.65)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -127L))

Plotting this signal over time gives me the following plot:
ggplot(df, aes(x = timestamp_zero, y = diff_YR)) +
  geom_line()

I want to do is select the values (i.e. diff_YR, timepoint, steering wheel angle (SWA) etc) at the peak of the diff_YR signal. At first I used the following code:
upperthreshold = 0.13 # upper threshold for consistent steering response
lowerthreshold = 0.05 # lower threshold for when response initiated

steering_peak <- df %>%
  dplyr::select(ppid_trialn, diff_YR, timestamp_zero, SWA) %>%
  group_by(ppid_trialn) %>%
  dplyr::mutate(peakYaw = max(diff_YR)) %>%
  dplyr::filter(max(diff_YR) > upperthreshold, min(diff_YR) < lowerthreshold) %>%
  slice(1:max(which(diff_YR > upperthreshold, 1))) %>%
  slice(max(which(diff_YR == peakYaw))) %>%
  ungroup() %>%
  transmute(ppid_trialn, PeakYawRateChange = diff_YR, PeakSteeringTime = timestamp_zero, peakSWA = SWA)

The above code groups by trial. It then identifies the maximum diff_YR signal. I  then filter the signal to include the largest diff_YR value over the upper threshold, and the smallest diff_YR below the lower threshold. Then I slice the diff_YR vector from the the first value to the max diff_YR value that is above the upper threshold. I then slice the row of the data when the diff_YR is equal to the maximum yaw.
This works fine if my maximum diff_YR value is at the peak however, as you can see from the plot, this is not the case for this trial. Hence I end up selecting a "peak" which is not the true peak of the signal. One suggestion I had was to filter the timecourse where the peak occurs. However for some trials, the peak is later in the timecourse thus I need all of it when looking at multiple trials.
Hence my question is: does anybody have an idea for how I could define and filter just the bellshapred curve response for each of my trials, and then select the row of data in my dataframe that corresponds to my true peak of the diff_YR signal.
Many thanks in advance! 

A:

You can try using packages that find peaks and allow you to define threshold etc, for example below, I use findpeaks from pracma, where you can provide a few options such as minimum peak height and minimum peak distance. I iterate through a few settings for minpeakpeakdistance, because it's easier to set something for minimum peak height:
library(pracma)
plotPeaks = function(df,i){
p = findpeaks(df$diff_YR,minpeakdistance=i,minpeakheight=0)
plot(df$timestamp_zero,df$diff_YR,main=paste("minpeakdistance=",i))
abline(v=df$timestamp_zero[p[,2]])
}

par(mfrow=c(2,2))
for(I in c(3,5,7,9)){
plotPeaks(df,I)
}

Now if you combine a minimum peak distance of 7 with a minimum peak height of say 0.1 you might get something sensible.
The find peak above does not take into account time, and it is also susceptible to bumps as you can see. You can also smooth the data first, or take a moving average to even out these:
df$mean_diff_YR = rollmean(df$diff_YR,5,fill=NA)
p = findpeaks(df$mean_diff_YR,minpeakdistance=7,minpeakheight=0)
plot(df$timestamp_zero,df$diff_YR,main=paste("minpeakdistance=",7))
abline(v=df$timestamp_zero[p[,2]])

