How to perform a polynomial equation and curve?

Question

I have this dataframe:

            id power     hr    fr    VE     VO2    VCO2  PETCO2 percent_VO2 percent_power
1  BM06-PRD-S1    25 119.25 18.25 19.00 0.61675 0.58225 37.6425    48.87084      25.00000
2  BM06-PRD-S1    40 126.00 18.00 20.75 0.71700 0.65950 39.2175    56.81458      40.00000
3  BM06-PRD-S1    55 133.50 20.75 25.00 0.86275 0.82750 41.2150    68.36371      55.00000
4  BM06-PRD-S1    70 147.25 18.25 29.00 0.98575 1.04550 41.7050    78.11014      70.00000
5  BM06-PRD-S1    85 158.50 22.25 39.25 1.13000 1.30525 41.1425    89.54041      85.00000
6  BM06-PRD-S1   100 168.75 27.75 51.00 1.26200 1.61150 38.8925   100.00000     100.00000
7  CB19-PRD-S1    25  98.75 18.50 25.00 0.88350 0.80475 40.7550    36.15715      13.15789
8  CB19-PRD-S1    40  98.25 20.00 25.50 0.94575 0.82900 41.4675    38.70473      21.05263
9  CB19-PRD-S1    55 102.00 19.75 28.50 1.08125 0.95800 42.2775    44.25005      28.94737
10 CB19-PRD-S1    70 107.50 20.50 34.25 1.24400 1.14275 42.6450    50.91058      36.84211
11 CB19-PRD-S1    85 111.00 21.25 35.50 1.30475 1.19925 43.3600    53.39677      44.73684
12 CB19-PRD-S1   100 117.25 21.50 40.25 1.47350 1.42225 44.2650    60.30284      52.63158
13 CB19-PRD-S1   115 123.00 22.75 47.00 1.67900 1.68475 44.6400    68.71291      60.52632
14 CB19-PRD-S1   130 129.50 24.50 52.50 1.79075 1.87950 44.3425    73.28627      68.42105
15 CB19-PRD-S1   145 135.50 25.25 59.50 1.96000 2.13525 44.7300    80.21281      76.31579
16 CB19-PRD-S1   160 145.25 26.75 64.50 2.04050 2.28350 43.8825    83.50726      84.21053
17 CB19-PRD-S1   175 151.25 30.50 83.00 2.34425 2.76050 41.6025    95.93820      92.10526
18 CB19-PRD-S1   190 161.75 33.75 92.25 2.44350 2.96850 40.0400   100.00000     100.00000
19 CC14-PRD-S1    20 102.50 19.00 18.25 0.59250 0.54825 37.7175    49.26211      22.22222
20 CC14-PRD-S1    30 110.25 18.75 19.75 0.66100 0.60325 38.5800    54.95739      33.33333
21 CC14-PRD-S1    40 113.25 18.50 20.75 0.74350 0.66025 39.2950    61.81667      44.44444
22 CC14-PRD-S1    50 122.50 20.00 23.50 0.87875 0.77325 40.5650    73.06173      55.55556
23 CC14-PRD-S1    60 126.25 17.50 26.25 0.94350 0.89375 41.3525    78.44523      66.66667
24 CC14-PRD-S1    70 132.00 16.50 28.00 0.99675 0.98525 42.7575    82.87258      77.77778
25 CC14-PRD-S1    80 145.00 18.50 32.75 1.11425 1.16275 42.5025    92.64186      88.88889
26 CC14-PRD-S1    90 153.50 19.50 37.25 1.20275 1.32700 42.0975   100.00000     100.00000
27 DA24-PRD-S1    25  88.00 18.50 15.75 0.53500 0.45075 37.2200    40.33170      21.73913
28 DA24-PRD-S1    40  93.25 18.50 16.25 0.58450 0.47775 38.3375    44.06332      34.78261
29 DA24-PRD-S1    55 103.75 19.00 20.25 0.76875 0.65450 40.1875    57.95326      47.82609
30 DA24-PRD-S1    70 119.00 20.75 28.00 0.98200 0.95525 41.5175    74.02940      60.86957
31 DA24-PRD-S1    85 133.25 22.75 34.75 1.09975 1.18325 41.4125    82.90614      73.91304
32 DA24-PRD-S1   100 145.00 27.50 45.75 1.25900 1.49700 39.1475    94.91142      86.95652
33 DA24-PRD-S1   115 155.25 36.50 64.75 1.32650 1.72500 33.0275   100.00000     100.00000

I am running a plot using ggplot and ggscatter:

ggplot(dftest, aes(percent_power, PETCO2)) +
  geom_point()

ggscatter(dftest, x = "percent_power", y = "PETCO2", add = "reg.line") +
  stat_cor(label.x = 20, label.y = 3.8) +
  stat_regline_equation(label.x = 20, label.y = 0.5) +
  xlab("Percentage of power (%)") + 
  geom_smooth(method = "lm", colour = "red") +
  ylab(expression(paste("PETC", O[2]," (mmHg)")))

I would like to perform a polynomial equation and curve because I am just able to run a linear regression.

Thank you!

Answer 1

I think what you're asking is how to how to fit a polynomial regression line to your data while still using the ggpubr functions to annotate it.

This is possible, but it seems that the in-built regression line can only be either a straight line or a loess model, neither of which is appropriate. However, you can fit a polynomial curve and get the equation and adjusted R-squared on the plot using the method below. In your case I have used a cubic formula, but you should choose your polynomial based on a known model or whatever makes most sense based on what you already know about the relationship between your variables.

You can use ggplot to add the actual line as suggested by @Roland, as long as this uses the same formula as the one you supply to stat_regline_equation

ggscatter(dftest, x = "percent_power", y = "PETCO2") +
  stat_regline_equation(label.x = 20, label.y = 0.5, 
                        formula = y ~ poly(x, 3),
                        aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")),) +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3)) +
  xlab("Percentage of power (%)") + 
  ylab(expression(paste("PETC", O[2]," (mmHg)")))

Which gives this result: