Search code examples
rplotbatch-processingnonlinear-functions

Plot quadratic regression with equation displayed


I am trying to create a pdf page conatining 6 plots (3 rows and 2 columns) using a for loop. I am able to create the plots but i cant seem to automate adding a regression line to each plot.

I am trying the following code.

#Dummy data
Data1 <- data.frame(flow = c(8,8.5,6,7.1,9), SP_elev = c(20,11,5,25,50))
Data2 <- data.frame(flow = c(7,7.2,6.5,8.2,8.5), SP_elev = c(13,15,18,25,19))
Data3 <- data.frame(flow = c(2,3,5,7,9), SP_elev = c(20,25,28,30,35))
Data4 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
Data5 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
Data6 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(22,23,25,27,29))

#Create Vector list 
dataframes = list("Data1" = Data1, 
                  "Data2" = Data2, 
                  "Data3" = Data3,
                  "Data4" = Data4,
                  "Data5" = Data5,
                  "Data6" = Data6) # I gave up here

# open the PDF device
pdf(file="Dummy_Example.pdf", paper="letter", height=10, width=8)

#Create array of plots 
par(mfrow=c(3,2))

#plot a with regression model
for (i in dataframes) {

plot (i[,c('flow', 'SP_elev')], xlab=expression(paste("Discharge (", ft^3, "/s)",sep = "")), ylab= "Elevation (m)", tck=0.02, adj = 0.5)

#plot regression curve
fit2<-lm(i$SP_elev ~ i$flow + I(i$flow^2), data=i) 
pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1] 
curve(pol2, lwd=1, add=T)
}

# close the PDF device
dev.off()

I got the regression curve code to work when i am individually producing a plot but it wont seem to work when i try to automate it.

In addition I also want to plot the equation of the regression curve on the graph.

Thanks,

dubbdan


Solution

  • I edited the code above to include the equation in the plot. Update: Now prettier equation.

    #Dummy data
    Data1 <- data.frame(flow = c(8,8.5,6,7.1,9), SP_elev = c(20,11,5,25,50))
    Data2 <- data.frame(flow = c(7,7.2,6.5,8.2,8.5), SP_elev = c(13,15,18,25,19))
    Data3 <- data.frame(flow = c(2,3,5,7,9), SP_elev = c(20,25,28,30,35))
    Data4 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
    Data5 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(13,15,18,25,19))
    Data6 <- data.frame(flow = c(1,4,6,8,9), SP_elev = c(22,23,25,27,29))
    
    #Create Vector list 
    dataframes = list("Data1" = Data1, 
                      "Data2" = Data2, 
                      "Data3" = Data3,
                      "Data4" = Data4,
                      "Data5" = Data5,
                      "Data6" = Data6) # I gave up here
    
    # open the PDF device
    pdf(file="Dummy_Example.pdf", paper="letter", height=10, width=8)
    
    #Create array of plots 
    par(mfrow=c(3,2))
    
    #plot a with regression model
    for (i in dataframes) {
    
    plot (SP_elev ~ flow, data=i,
          xlab=expression(paste("Discharge (", ft^3, "/s)",sep = "")),
          ylab= "Elevation (m)", tck=0.02, adj = 0.5)
    
    #plot regression curve
    fit2<-lm(SP_elev ~ flow + I(flow^2), data=i) 
    pol2 <- function(x) fit2$coefficient[3]*x^2 + fit2$coefficient[2]*x + fit2$coefficient[1] 
    curve(pol2, lwd=1, add=T, col="blue")
    
    xm <- min(i$flow)
    ym <- max(i$SP_elev)
    
    co <- signif(coef(fit2),3)
    text(xm, ym,
     bquote(y==.(co[3])*x^2 + .(co[2])*x + .(co[1])),
     adj=c(0,1))
    
    
    }
    dev.off()