Adding a trendline in basic R code (no GGPlot) for multiple plot graph of exponential and 2 factor polynomial curves

by

I am trying to add trendlines to a multi-plot graph using basic R code. Could use some help to to add trendlines for non-linear functions, and am not sure if this is possible with basic R. I may need to upgrade my skillset to ggplot and thought I would check here first.

I have soil data for two locations, a reference site, and a site grazed by livestock. There are four litter decomposition metrics that I’ve tracked over time at each location: S, k, Bg and Br. The purpose of this experiment was to determine the shape of the decomposition curve for each metric. R2 analysis identified that S has best fit to an exponential model of decay, while k, Bg and Br are best fit to a 2 factor polynomial curve.

I’m trying to visualize the shapes of each best fit curve to each metric over time. I have tried to add trendlines but can’t make anything work that isn’t linear. I have tried different abline functions but am not sure this is the right approach. Thank you Allan Cameron for your response. I’ve edited this question to provide the full code for each of my models. I’m not sure how to upload my dataframe but have also provided some header data from my .csv.

This is as far as I can get visualizing the four plots

Here is the R code that I’ve been working with so far:

# Load the farm data file
soil_data <- read.csv("C1Sandown020223.csv", header = TRUE)

# Check the header of my data
head(soil_data)

# Look at the whole dataset
soil_data

# Create a new plot
plot.new()

# Set up multiple plots to visualize data (4 rows, 2 columns)
par(mfrow = c(1, 1))

# Plot S for each site

# Graph linear model of S for Reference against the number of 
# days in the ground
plot(soil_data$Days1, soil_data$S1, col = "darkgreen", xlab = "Days",
    ylab = "S", ylim = c(0, 0.7), las = 1,
   main = "Reference")

modelS1 <- lm(S1 ~ Days1, data = soil_data)
abline(modelS, col = "green")
legend(100, 0.70, bty="n", legend=c("R2 = 0.873"))
summary(modelS1)

# Graph linear model of S for Grazed against the number of 
# days in the ground
plot(soil_data$Days2, soil_data$S2, col = "darkgreen", xlab = "Days",
    ylab = "S", ylim = c(0, 0.7), las = 1,
     main = "Grazed")

modelS2 <- lm(S2 ~ Days2, data = soil_data)
abline(modelS2, col = "green")
legend(100, 0.70, bty="n", legend=c("R2 = 0.839"))
summary(modelS2)

# Graph linear model of S for Combined Sites against the 
# number of days in the ground
plot(soil_data$Days, soil_data$S, col = "darkblue", xlab = "Days",
    ylab = "S", ylim = c(0, 0.7), las = 1,
 main = "Combined")

modelS <- lm(S ~ Days, data = soil_data)
abline(modelS, col = "darkblue")
legend(100, 0.70, bty="n", legend=c("R2 = 0.857"))
summary(modelS)

# Plot linear model of k for each site 

# Graph linear model of k for Reference against the number of 
# days in the ground
plot(soil_data$Days1, soil_data$k1, col = "darkred", xlab = "Days",
 ylab = "k", ylim = c(0, 0.05), las = 1,
 main = "Reference")

modelk1 <- lm(k1 ~ Days1, data = soil_data)
abline(modelk1, col = "darkred")
legend("topright", bty="n", legend=c("R2 = 0.523"))
summary(modelk1)

# Graph linear model of k for Grazed against the # days in the ground
plot(soil_data$Days2, soil_data$k2, col = "darkred", xlab = "Days",
 ylab = "k", ylim = c(0, 0.05), las = 1,
 main = "Grazed")

modelk2 <- lm(k2 ~ Days2, data = soil_data)
abline(modelk2, col = "darkred")
legend("topright", bty="n", legend=c("R2 = 0.459"))
summary(modelk2)

# Graph linear model of k for Combined Sites against the # days 
plot(soil_data$Days, soil_data$k, col = "darkblue", xlab = "Days",
 ylab = "k", ylim = c(0, 0.05), las = 1,
 main = "Combined")

modelk <- lm(k ~ Days, data = soil_data)
abline(modelk, col = "darkblue")
legend("topright", bty="n", legend=c("R2 = 0.465"))
summary(modelk)

# Plot linear model of Bg for each site 

# Graph linear model of Bg for Reference against the number of days 
plot(soil_data$Days1, soil_data$BG1, col = "green", xlab = "Days",
 ylab = "Bg", ylim = c(0.3, 0.8), las = 1,
 main = "Reference")

modelBG1 <- lm(BG1 ~ Days1, data = soil_data)
abline(modelBG1, col = "green")
legend(100, 0.40, bty="n", legend=c("R2 = 0.873"))
summary(modelBG1)

# Graph linear model of Bg for Grazed against the number of days 
plot(soil_data$Days2, soil_data$BG2, col = "green", xlab = "Days",
 ylab = "Bg", ylim = c(0.3, 0.8), las = 1,
 main = "Grazed")

modelBG2 <- lm(BG2 ~ Days2, data = soil_data)
abline(modelBG2, col = "green")
legend(100, 0.40, bty="n", legend=c("R2 = 0.839"))
summary(modelBG2)

# Graph linear model of Bg for Combined Sites against # of days 
plot(soil_data$Days, soil_data$BG, col = "darkblue", xlab = "Days",
 ylab = "Bg", ylim = c(0.3, 0.8), las = 1,
 main = "Combined")

modelBG <- lm(BG ~ Days, data = soil_data)
abline(modelBG, col = "darkblue")
legend(100, 0.40, bty="n", legend=c("R2 = 0.857"))
summary(modelBG)

# Plot linear model of Br for each site

# Graph linear model of Br for reference against the # days 
plot(soil_data$Days1, soil_data$BR1, col = "darkred", xlab = "Days",
 ylab = "Br", ylim = c(0.2, 0.5), las = 1,
 main = "Reference")

modelBR1 <- lm(BR1 ~ Days1, data = soil_data)
abline(modelBR1, col = "darkred")
legend("bottomright", bty="n", legend=c("R2 = 0.921"))
summary(modelBR1)

# Graph linear model of Br for Grazed against the # of days 
plot(soil_data$Days2, soil_data$BR2, col = "darkred", xlab = "Days",
 ylab = "Br", ylim = c(0.2, 0.5), las = 1,
 main = "Grazed")

modelBR2 <- lm(BR2 ~ Days2, data = soil_data)
abline(modelBR2, col = "darkred")
legend("bottomright", bty="n", legend=c("R2 = 0.833"))
summary(modelBR2)

# Graph linear model of Br for Combined against the number of days 
plot(soil_data$Days, soil_data$BR, col = "darkblue", xlab = "Days",
 ylab = "Br", ylim = c(0.2, 0.5), las = 1,
 main = "Combined")

modelBR <- lm(BR ~ Days, data = soil_data)
abline(modelBR, col = "darkblue")
legend("bottomright", bty="n", legend=c("R2 = 0.885"))
summary(modelBR)

# EXPONENTIAL MODELS FOR EACH MATERIAL AND MEASURE 
# (log value of y axis)
# take the log of response variable, plot versus predictor 
# variable (days)
# method eg. https://www.statology.org/exponential-regression-in-r/

# Exponential model of S

# S Reference Exp
expSRef = lm(log(soil_data$S1) ~ soil_data$Days1)
summary(expSRef)

# S Grazed Exp
expSGraz = lm(log(soil_data$S2) ~ soil_data$Days2)
summary(expSGraz)

# S Combined Exp
expS = lm(log(soil_data$S) ~ soil_data$Days)
summary(expS)

# Exponential model of k for each site 

# k Reference Exp
expkRef = lm(log(soil_data$k1) ~ soil_data$Days1)
summary(expkRef)

# k Grazed Exp
expkGraz = lm(log(soil_data$k2) ~ soil_data$Days2)
summary(expkGraz)

# k Combined Exp
expk = lm(log(soil_data$k) ~ soil_data$Days)
summary(expk)

# Exponential model of Bg for each site 

# Bg Reference Exp
expBgRef = lm(log(soil_data$BG1) ~ soil_data$Days1)
summary(expBgRef)

# Bg Grazed Exp
expBgGraz = lm(log(soil_data$BG2) ~ soil_data$Days2)
summary(expBgGraz)

# Bg Combined Exp
expBg = lm(log(soil_data$BG) ~ soil_data$Days)
summary(expBg)

# Exponential model of Br for each site

# Br Reference Exp
expBrRef = lm(log(soil_data$BR1) ~ soil_data$Days1)
summary(expBrRef)

# Br Grazed Exp
expBrGraz = lm(log(soil_data$BR2) ~ soil_data$Days2)
summary(expBrGraz)

# Br Combined Exp
expBr = lm(log(soil_data$BR) ~ soil_data$Days)
summary(expBr)


# LOGARITHMIC MODELS FOR EACH MATERIAL AND MEASURE  
# (log value of x axis)
# https://www.statology.org/logarithmic-regression-in-r/

# Logarithmic model of S

# S Reference Log
logSRef = lm(log(soil_data$Days1) ~ soil_data$S1)
summary(logSRef)

# S Grazed Log
logSGraz = lm(log(soil_data$Days2) ~ soil_data$S2)
summary(logSGraz)

# S Combined Log
logS = lm(log(soil_data$Days) ~ soil_data$S)
summary(logS)


# Logarithmic model of k for each site 

# k Reference Log
logkRef = lm(log(soil_data$Days1) ~ soil_data$k1)
summary(logkRef)

# k Grazed Log
logkGraz = lm(log(soil_data$Days2) ~ soil_data$k2)
summary(logkGraz)

# k Combined Log
logk = lm(log(soil_data$Days) ~ soil_data$k)
summary(logk)

# Logarithmic model of Bg for each site 

# Bg Reference Log
logBgRef = lm(log(soil_data$Days1) ~ soil_data$BG1)
summary(logBgRef)

# Bg Grazed Log
logBgGraz = lm(log(soil_data$Days2) ~ soil_data$BG2)
summary(logBgGraz)

# Bg Combined Log
logBg = lm(log(soil_data$Days) ~ soil_data$BG)
summary(logBg)

# Logarithmic model of Br for each site

# Br Reference Log
logBrRef = lm(log(soil_data$Days1) ~ soil_data$BR1)
summary(logBrRef)

# Br Grazed Log
logBrGraz = lm(log(soil_data$Days2) ~ soil_data$BR2)
summary(logBrGraz)

# Br Combined Log
logBr = lm(log(soil_data$Days) ~ soil_data$BR)
summary(logBr)


# POLYNOMIAL MODEL OF EACH METRIC AND MEASURE - 2 dimensions

# Polynomial model of S

# S Reference Poly 2 dimensions
polySRef = lm(S1~poly(Days1, 2, raw = TRUE), data=soil_data)
summary(polySRef)

# S Grazed Poly 2 dimensions
polySGraz = lm(S2~poly(Days2, 2, raw = TRUE), data=soil_data)
summary(polySGraz)

# S Combined Poly 2 dimensions
polyS = lm(S~poly(Days, 2, raw = TRUE), data=soil_data)
summary(polyS)


# Polynomial model of k for each site 

# k Reference Poly 2 dimensions
polykRef = lm(k1~poly(Days1, 2, raw = TRUE), data=soil_data)
summary(polykRef)

# k Grazed Poly 2 dimensions
polykGraz = lm(k2~poly(Days2, 2, raw = TRUE), data=soil_data)
summary(polykGraz)

# k Combined Poly 2 dimensions
polyk = lm(k~poly(Days, 2, raw = TRUE), data=soil_data)
summary(polyk)

# Polynomial model of Bg for each site 

# Bg Reference Poly 2 dimensions
polyBgRef = lm(BG1~poly(Days1, 2, raw = TRUE), data=soil_data)
summary(polyBgRef)

# Bg Grazed Poly 2 dimensions
polyBgGraz = lm(BG2~poly(Days2, 2, raw = TRUE), data=soil_data)
summary(polyBgGraz)

# Bg Combined Poly 2 dimensions
polyBg = lm(BG~poly(Days, 2, raw = TRUE), data=soil_data)
summary(polyBg)

# Polynomial model of Br for each site

# Br Reference Poly 2 dimensions
polyBrRef = lm(BR1~poly(Days1, 2, raw = TRUE), data=soil_data)
summary(polyBrRef)

# Br Grazed Poly 2 dimensions
polyBrGraz = lm(BR2~poly(Days2, 2, raw = TRUE), data=soil_data)
summary(polyBrGraz)

# Br Combined Poly 2 dimensions
polyBr = lm(BR~poly(Days, 2, raw = TRUE), data=soil_data)
summary(polyBr)


# Summary R squared data for each model

# LINEAR
summary(modelS1)$r.squared
summary(modelS2)$r.squared
summary(modelS)$r.squared
summary(modelk1)$r.squared
summary(modelk2)$r.squared
summary(modelk)$r.squared
summary(modelBG1)$r.squared
summary(modelBG2)$r.squared
summary(modelBG)$r.squared
summary(modelBR1)$r.squared
summary(modelBR2)$r.squared
summary(modelBR)$r.squared

# EXPONENTIAL
summary(expSRef)$r.squared
summary(expSGraz)$r.squared
summary(expS)$r.squared
summary(expkRef)$r.squared
summary(expkGraz)$r.squared
summary(expk)$r.squared
summary(expBgRef)$r.squared
summary(expBgGraz)$r.squared
summary(expBg)$r.squared
summary(expBrRef)$r.squared
summary(expBrGraz)$r.squared
summary(expBr)$r.squared

# LOGARITHMIC 
summary(logSRef)$r.squared
summary(logSGraz)$r.squared
summary(logS)$r.squared
summary(logkRef)$r.squared
summary(logkGraz)$r.squared
summary(logk)$r.squared
summary(logBgRef)$r.squared
summary(logBgGraz)$r.squared
summary(logBg)$r.squared
summary(logBrRef)$r.squared
summary(logBrGraz)$r.squared
summary(logBr)$r.squared

# POLYNOMIAL - 2 dimensions
summary(polySRef)$r.squared
summary(polySGraz)$r.squared
summary(polyS)$r.squared
summary(polykRef)$r.squared
summary(polykGraz)$r.squared
summary(polyk)$r.squared
summary(polyBgRef)$r.squared
summary(polyBgGraz)$r.squared
summary(polyBg)$r.squared
summary(polyBrRef)$r.squared
summary(polyBrGraz)$r.squared
summary(polyBr)$r.squared

# Create combined dataset of S, k, Br, Bg, with 
# different shapes for each site

# Create a new plot
plot.new()

# Set up multiple plots to visualize data (4 rows, 2 columns)
par(mfrow = c(2, 2))

# Graph linear model of S for Combined Sites against # days 

plot(soil_data$Days[soil_data$Type == "Grazed"], 
soil_data$S[soil_data$Type == "Grazed"], 
 pch = 16, cex = 1, col = "darkblue", ylim = c(0, 0.7), 
xlab = "Days", ylab = "S", main = "S")
points(soil_data$Days[soil_data$Type == "Reference"], 
soil_data$S[soil_data$Type == "Reference"], 
   pch = 17, cex = 1, col = "darkred")
# this abline function is giving a straight line and 
# not the actual model, does the exp not make 
# an exponential trendline? 
abline(expSGraz, col = "darkblue")
abline(expSRef, col = "darkred")
legend("topright", legend = unique(soil_data$Type), pch = c(16, 17),
   col = c("darkblue", "darkred"), cex = 0.8, title = "Data Type")


# Now for k

plot(soil_data$Days[soil_data$Type == "Grazed"], 
soil_data$k[soil_data$Type == "Grazed"], 
 pch = 16, cex = 1, col = "darkblue", xlab = "Days", 
ylab = "k", main = "k")
points(soil_data$Days[soil_data$Type == "Reference"],
soil_data$k[soil_data$Type == "Reference"], 
   pch = 17, cex = 1, col = "darkred")
# this abline function is giving a straight line and 
# not the actual model 
abline(polykGraz, col = "darkblue")
abline(polykRef, col = "darkred")
legend("topright", legend = unique(soil_data$Type), pch = c(16, 17),
   col = c("darkblue", "darkred"), cex = 0.8, title = "Data Type")

# Bg

plot(soil_data$Days[soil_data$Type == "Grazed"],
soil_data$BG[soil_data$Type == "Grazed"], 
 pch = 16, cex = 1, col = "darkblue", xlab = "Days", 
ylab = "Bg", main = "Bg")
points(soil_data$Days[soil_data$Type == "Reference"],
soil_data$BG[soil_data$Type == "Reference"], 
   pch = 17, cex = 1, col = "darkred")
# this abline function is giving a straight line and 
# not the actual model 
abline(polyBgGraz, col = "darkblue")
abline(polyBgRef, col = "darkred")
legend("topleft", legend = unique(soil_data$Type), pch = c(16, 17),
   col = c("darkblue", "darkred"), cex = 0.8, title = "Data Type")

#Br

plot(soil_data$Days[soil_data$Type == "Grazed"],
soil_data$BR[soil_data$Type == "Grazed"], 
 pch = 16, cex = 1, col = "darkblue", xlab = "Days", 
ylab = "Br", main = "Br")
points(soil_data$Days[soil_data$Type == "Reference"],
soil_data$BR[soil_data$Type == "Reference"], 
   pch = 17, cex = 1, col = "darkred")
# this abline function is giving a straight line 
# and not the actual model 
abline(polyBrGraz, col = "darkblue")
abline(polyBrRef, col = "darkred")

legend("topleft", legend = unique(soil_data$Type), pch = c(16, 17),
   col = c("darkblue", "darkred"), cex = 0.8, title = "Data Type")

My Data

dput(soil_data)
structure(list(Site = c(23L, 19L, 18L, 16L, 20L, 17L, 21L, 22L,
24L, 6L, 7L, 2L, 5L, 3L, 12L, 13L, 8L, 14L, 9L, 15L, 4L, 1L,
10L, 11L), Type = c(“Grazed”, “Grazed”, “Grazed”, “Grazed”, “Grazed”,
“Grazed”, “Grazed”, “Grazed”, “Grazed”, “Reference”, “Reference”,
“Reference”, “Reference”, “Reference”, “Reference”, “Reference”,
“Reference”, “Reference”, “Reference”, “Reference”, “Reference”,
“Reference”, “Reference”, “Reference”), Days = c(35L, 49L, 63L,
77L, 91L, 98L, 105L, 119L, 133L, 35L, 42L, 49L, 56L, 63L, 70L,
77L, 84L, 91L, 98L, 105L, 112L, 119L, 126L, 133L), S = c(0.566,
0.593, 0.57, 0.285, 0.225, 0.197, 0.269, 0.205, 0.119, 0.573,
0.566, 0.42, 0.628, 0.507, 0.455, 0.268, 0.306, 0.231, 0.196,
0.198, 0.164, 0.101, 0.127, 0.114), k = c(0.0591, 0.0089, 0.0281,
0.0113, 0.0061, 0.0099, 0.0052, 0.0106, 0.0079, 0.0484, 0.0227,
0.0216, 0.0116, 0.0237, 0.0186, 0.0136, 0.0227, 0.01, 0.0105,
0.0083, 0.0082, 0.0125, 0.0124, 0.0072), BG = c(0.365, 0.343,
0.362, 0.602, 0.653, 0.676, 0.615, 0.669, 0.742, 0.359, 0.365,
0.488, 0.313, 0.415, 0.459, 0.616, 0.584, 0.647, 0.677, 0.676,
0.704, 0.757, 0.735, 0.746), BR = c(0.239, 0.225, 0.366, 0.395,
0.428, 0.443, 0.403, 0.439, 0.486, 0.235, 0.239, 0.32, 0.325,
0.327, 0.301, 0.404, 0.383, 0.433, 0.444, 0.443, 0.461, 0.496,
0.482, 0.489), Days2 = c(35L, 49L, 63L, 77L, 91L, 98L, 105L,
119L, 133L, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA), S2 = c(0.566, 0.593, 0.57, 0.285, 0.225, 0.197, 0.269,
0.205, 0.119, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA), k2 = c(0.0591, 0.0089, 0.0281, 0.0113, 0.0061, 0.0099,
0.0052, 0.0106, 0.0079, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA), BG2 = c(0.365, 0.343, 0.362, 0.602, 0.653,
0.676, 0.615, 0.669, 0.742, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA), BR2 = c(0.239, 0.225, 0.366, 0.395,
0.428, 0.443, 0.403, 0.439, 0.486, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA), Days1 = c(NA, NA, NA, NA, NA,
NA, NA, NA, NA, 35L, 42L, 49L, 56L, 63L, 70L, 77L, 84L, 91L,
98L, 105L, 112L, 119L, 126L, 133L), S1 = c(NA, NA, NA, NA, NA,
NA, NA, NA, NA, 0.573, 0.566, 0.42, 0.628, 0.507, 0.455, 0.268,
0.306, 0.231, 0.196, 0.198, 0.164, 0.101, 0.127, 0.114), k1 = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, 0.0484, 0.0227, 0.0216, 0.0116,
0.0237, 0.0186, 0.0136, 0.0227, 0.01, 0.0105, 0.0083, 0.0082,
0.0125, 0.0124, 0.0072), BG1 = c(NA, NA, NA, NA, NA, NA, NA,
NA, NA, 0.359, 0.365, 0.488, 0.313, 0.415, 0.459, 0.616, 0.584,
0.647, 0.677, 0.676, 0.704, 0.757, 0.735, 0.746), BR1 = c(NA,
NA, NA, NA, NA, NA, NA, NA, NA, 0.235, 0.239, 0.32, 0.325, 0.327,
0.301, 0.404, 0.383, 0.433, 0.444, 0.443, 0.461, 0.496, 0.482,
0.489)), class = “data.frame”, row.names = c(NA, -24L))

  • 2

    Hi EBH. Welcome to Stack Overflow. It’s going to be difficult to give you an answer with the information you have provided. We don’t have soil_data, nor the code you used to create your models, so we cannot test or verify solutions. As a result, any code we offer could be wide of the mark and there is no way we can demonstrate that it works. Please make your question reproducible. If you edit your question to include the output from soil_data, and the code used to generate your models, it will be easy to help. Without that, we can only give vague guidance in the comments. Thanks.

    – 

  • Thank you @AllanCameron! I have updated the question. Please let me know if this is enough information to support a response. Clearly I’m new to this 😐

    – 

  • thank you. But can you please include the data as actual text that we can copy-paste into our consoles (perhaps the output from dput(soil_data)? Obviously, we can’t run your code with a picture of the data – we need the data itself (unless someone decides to transcribe your picture back into data, but this seems like an unreasonable thing to expect!)

    – 

  • Thank you for the suggested prompt to pull the data out into this form. I’ve added it to the question – let me know if this works

    – 

Leave a Comment