Obtain X value at which Max Y value occurs from spline model and extrapolate to find where model intersects x-axis

General
1

Hello

I am working with data generated by respiration rate in response to temperature.

I want to fit a smooth spline and natural spline to my data (I already did this), and from either and/or both of those models determine the X value (represented by X5) at which the MAX Y value (represented by X4) occurs. Once extracted I want to display the MAX point on my plot (NB I already know how to find the MAX Y and corresponding X using just the data.frame
I also want to extrapolate my data to find where the model intersects the x-axis on either the left or right side of the plot.

I gave my example data and, my code fitting the splines and obtaining MAX Y using the vector from my data below.

PS Ignore most of those packages

library(zoo)
#> 
#> Attaching package: 'zoo'
#> The following objects are masked from 'package:base':
#> 
#>     as.Date, as.Date.numeric
library(lubridate)
#> 
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#> 
#>     date, intersect, setdiff, union
library(tidyverse)
library(readxl)
library(readr)
library(ggplot2)
library(pacman)
library(reshape2)
#> 
#> Attaching package: 'reshape2'
#> The following object is masked from 'package:tidyr':
#> 
#>     smiths
library(splines)

##############################################################################################################


# create dataframe
X5 <- c(20.0857313,
       20.2006048,
       20.33861734,
       20.48452669,
       20.63168701,
       20.79547963,
       20.98120269,
       21.16778097,
       21.36226786,
       21.56228305,
       21.78009325,
       21.9907538,
       22.21212998,
       22.43460528,
       22.64577427,
       22.87884339,
       23.11815104,
       23.35915888,
       23.60532229,
       23.84629768,
       24.08781313,
       24.34224763,
       24.59094248,
       24.84459757,
       25.07997528,
       25.31870461,
       25.56135172,
       25.79729375,
       26.03642188,
       26.27510938,
       26.52024513,
       26.76156151,
       27.00034359,
       27.2462221,
       27.51208148,
       27.75439341,
       27.99372057,
       28.23603718,
       28.46266882,
       28.67969867,
       28.87219326,
       29.07419544,
       29.29385518,
       29.53276831,
       29.77350961,
       30.0348893,
       30.30289255,
       30.53915662,
       30.80427146,
       31.09287906,
       31.41811948,
       31.74206118,
       32.0393169,
       32.3057906,
       32.58206629,
       32.87828599,
       33.16526798,
       33.42830529,
       33.71557346,
       33.9922445,
       34.23972463,
       34.49415729,
       34.7722132,
       35.02893689,
       35.28454039,
       35.5389079,
       35.8231931,
       36.10830968,
       36.39164222,
       36.65346271,
       36.90514364,
       37.20379735,
       37.46584994,
       37.73181209,
       38.01365826,
       38.31509393,
       38.60053864,
       38.91812969,
       39.22217979,
       39.50730872,
       39.802358,
       40.08579487,
       40.33521347,
       40.57716969,
       40.79089089,
       41.04978681,
       41.34364629,
       41.62344038,
       41.86821976,
       42.16751008,
       42.46597745,
       42.72254633,
       43.01346491,
       43.3221519,
       43.60901729,
       43.83435099,
       44.10540511,
       44.39115612,
       44.68864736,
       44.98294805,
       45.29779589,
       45.57835765,
       45.89134984,
       46.18074239,
       46.45019306,
       46.7357955,
       47.05500415,
       47.32244008,
       47.64222648,
       47.99330107,
       48.3117349,
       48.60130314,
       48.91275987,
       49.26994847)

X4 <- c(44.07028337,
44.30479187,
41.84341936,
45.45282442,
44.7268717,
44.77605349,
48.07465638,
46.54220775,
45.27858738,
45.65923262,
43.35510699,
48.67790975,
49.63068536,
49.1343319,
50.08969849,
50.5089601,
49.50857752,
53.44111837,
55.38744667,
56.62011579,
55.81142272,
56.44124049,
60.29871937,
61.42926584,
60.07091676,
59.89739533,
64.30068066,
66.34449833,
67.57523983,
67.32368324,
71.64170341,
71.51606386,
70.84460881,
70.25193327,
74.60831589,
73.11810445,
72.65562864,
73.81531117,
74.48278155,
75.98856074,
73.22999998,
74.66923839,
74.3257666,
81.36855953,
80.54377644,
82.45988801,
83.30434199,
76.23624,
87.57180724,
88.16171033,
88.58291481,
90.1037515,
88.69838003,
88.79216561,
92.21198308,
92.41351747,
91.46948895,
89.77822239,
89.06744027,
88.47788658,
86.0380284,
88.32497018,
89.80440745,
92.79966501,
93.70495761,
89.61016982,
93.81707579,
88.9454452,
93.31642691,
91.40727389,
92.22236809,
93.96125291,
93.70331553,
91.52366664,
91.58078309,
92.19943323,
92.5447887,
92.31490548,
90.40022368,
93.0475149,
92.91598461,
90.19452399,
89.45602688,
89.84914851,
88.61844803,
89.16295867,
89.93636587,
89.08043077,
87.25242338,
87.20588203,
87.89543603,
87.97789447,
85.98634579,
85.36283357,
86.56698373,
87.49965799,
88.65618018,
89.30224619,
88.44160182,
92.59997657,
92.59414076,
86.83093222,
81.00150256,
66.62084299,
60.67011001,
64.5729937,
57.54471422,
58.69267483,
56.29551649,
52.43103495,
45.95466922,
38.58173711,
33.81493804,
32.27193129)

IMP2 <- data.frame(X5, X4)

#alternative object to store transformed data
df2 <- IMP2

#Transform data
df2$X4 <- log(df2$X4)
df2$X5 <- 1/(df2$X5+273.15)*1000 

# fitting splines
#smooth splines
SS2 <- smooth.spline(df2$X5, df2$X4, cv = TRUE)
plot(df2$X5, df2$X4)
lines(SS2, col=2, lwd=3)



#natural splines
NS2 <- lm(df2$X4 ~ ns(df2$X5, df = 6))
plot(df2$X5, df2$X4)
lines(df2$X5, predict(NS2), col=2, lwd=3)


#find max from vectors, not using spline model
which.max(df2$X4)
#> [1] 72
df2$X5[which.max(df2$X4)]
#> [1] 3.222129
abline(v=df2$X5[which.max(df2$X4)], col="red")

Created on 2021-05-04 by the reprex package (v1.0.0)

2

Welcome to the community. For natural spline, I can find the X5 for which prediction is max as below. Please let me know if I understand correctly.

lmax<-df2$X5[which.max(predict(NS2))]
abline(v=lmax, col="blue")
1 Like
3

@mhakanda Yes that works perfectly, thank you.
I also want to extrapolate the natural spline (NS2) model to find where the model intersects the x-axis (determine the x-value here) on either the left or right side of the plot

4

Getting roots is kind of guess the interval containing roots and try it out. We see there are two roots, in left and right. First we plot an extended plot, guess it and then confirm it.

# For natural spline, it is the same model as your one
NS3<- lm(X4 ~ ns(X5, df = 6), data = df2)
#choose the interval and increment
s1<-seq(2.8, 4.2, 0.001)
p1<-predict(NS3, newdata = data.frame(X5 = s1))

# see the plot and guess where roots might be
plot(p1~s1, col=2, lwd=3, type = "l", xlab = "X5", ylab = "X4")
grid(NA, 5, lwd = 2)
# so it seems there are two zeors, 
# one around 3 and other one around 4

# so, to get first one, choose interval 1.8 to 3.4
s11<-seq(2.8,3.4, .001)
r1<-s11[which.min(abs(predict(NS3, newdata = data.frame(X5 = s11))))]
r1
#> [1] 2.929
# so, to get second one, choose interval 3.4to 4.2
s22<-seq(3.4,4.2, .001)
r2<-s22[which.min(abs(predict(NS3, newdata = data.frame(X5 = s22))))]
r2
#> [1] 3.992

abline(v=r1, col="blue")
abline(v=r2, col="blue")

1 Like
5

@mhakanda Perfect, thank you so much for all your help.

6

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.