polynomial.plot.R

#title: code sample for ploting polynomial model  

# load packages
library(readxl) # load readxl package
library(ggplot2) # load ggplot2 paackage
library(dplyr) # load dplyr package
library(polynom) # load polynom
library(gsheet) # load gsheet 

# load data 
# x is nitrogen and y is yield
data <- gsheet2tbl("docs.google.com/spreadsheets/d/10v0S4fUh0zylHSVVA7mX2bEIjNfm2Vs49tpwXm0fFBg")

# plot scatter 
p <- data %>% ggplot(aes(x = x, y = y)) +  geom_point(alpha = 2/10, shape = 21, 
                                                      fill = "blue", color = "black", size = 5)
p # showing graph

# check models
tidy(lm(y ~ poly(x, 2, raw=TRUE), data = data)) # polynomial power2
tidy(lm(y ~ poly(x, 3, raw=TRUE), data = data)) # polynomial power3

# plot graph with trand line
p1 <- p + stat_smooth(method="lm", se=TRUE, fill=NA, 
                      formula = y ~ poly(x, 2, raw=TRUE),colour="red")

# extract the value from model
my.formula <- y ~ poly(x, 2, raw = TRUE)

m <- lm(my.formula, data)

my.eq <- as.character(signif(as.polynomial(coef(m)), 2))

label.text <- paste(gsub("x", "~italic(x)", my.eq, fixed = TRUE),
                    paste("italic(R)^2",  
                          format(summary(m)$r.squared, digits = 2), 
                          sep = "~`=`~"),
                    sep = "~~~~")

# plot graph with equation
p1 + annotate("text", x = 10, y = 750, label = label.text, hjust=0, size=8, 
             family="Times", parse=TRUE)

# save p-value report 
report <- tidy(lm(y ~ poly(x, 2, raw=TRUE), data = data))

write.csv(file = "report.csv", report)