MLRpack: Introductory Pacakge

# Load libraries needed
library(shiny)
#library(ggplot2)
library(purrr)
library(rootSolve)
library(scatterplot3d)
source("Rcode.r")
library(dplyr)

spruce.df = read.csv("SPRUCE.csv")

d = spruce.df$BHDiameter



# Define UI for application that draws a histogram
ui <- fluidPage(

   # Application title
   titlePanel("Spruce Data Set: Piecewise Regression"),

   # Sidebar with a slider input for number of bins

   sidebarLayout(
      sidebarPanel(
         sliderInput("xk",
                     "Choose first knot xk1.",
                     min = min(d),
                     max = max(d),
                     value = 17.44165,
                     step=0.01),

         sliderInput("xk2",
                     "Choose second knot xk2. Note this must be done after xk1 is chosen.",
                     min = min(d),
                     max = max(d),
                     value = 17.44165,
                     step=0.01),

         sliderInput("delta",
                     "Choose a value for delta. This controls the width of the intervals around
                     the knots. The intervals have length 2*Delta. Sliding changes the last table.",
                     min = 0.0,
                     max = 5.0,
                     value = c(1.0),
                     step=0.1)

      ),

      # Show a plot of the generated distribution
      mainPanel(
         plotOutput("regressPlot"),
         plotOutput("R2"),
         tableOutput("knots"),
         # table of data
         tableOutput("tab"),
         plotOutput("allroots")
      )
   )
)

# Define server logic required to draw a histogram
server <- function(input, output,session) {




   output$regressPlot <- renderPlot({
      plot(spruce.df,main="Piecewise regression",pch=21,bg="black")


      sp2.df=within(spruce.df, X<-(BHDiameter-input$xk)*(BHDiameter>input$xk))

      sp2.df=within(sp2.df, K<-(BHDiameter-input$xk2)*(BHDiameter>input$xk2))

      lmp = lm(Height ~ BHDiameter + X + K, data = sp2.df)
      tmp=summary(lmp) # tmp holds the summary info



      curve(myf2(x,xk=input$xk,xk2=input$xk2,coef=tmp$coefficients[,"Estimate"] ),
            add=TRUE,
            lwd=2,
            col="Blue")


      points(input$xk,myf(input$xk,input$xk,coef=tmp$coefficients[,"Estimate"] ),col="black",pch=21,bg="green",cex=2)

      points(input$xk2,myf2(input$xk2,input$xk,input$xk2,coef=tmp$coefficients[,"Estimate"] ),col="black",pch=21,bg="yellow",cex=2)

      #points(uroot()$root,myf(uroot()$root,uroot()$root,coef=tmp$coefficients[,"Estimate"] ),col="black",pch=21,bg="purple",cex=2)

      text(input$xk,16,
           paste("R sq.=",round(tmp$r.squared,4) ))

   })

   observe({
      val <- input$xk
       #Control the value, min, max, and step.
       #Step size is 2 when input value is even; 1 when value is odd.
      updateSliderInput(session, "xk2", value = val, step = .01)

   })

   #uroot = reactive({
    #  intv = input$intervalroot
     # uniroot(f=rsqdash, interval=intv, h=0.001,data=spruce.df, extendInt = "yes" )
   #})

   #urootall = reactive({
    #  intv = input$intervalroot
     # uniroot.all(f=rsqdash, interval=intv, h=0.001,data=spruce.df )
   #})

#   output$R2 <- renderPlot({
#      dsmooth = seq(min(d),max(d),length=1000)
#      r2=map_dbl(dsmooth, ~rsq(.x,data=spruce.df))
#      plot(dsmooth,r2,pch=21,bg="purple",
#           ylab = expression(R^2 ),
#           xlab = "knot",
#           cex=0.5,
#           main="Determination of possible knots",
#           type="p",
#           ylim = c(min(r2),1.1*max(r2)))
#      intv = input$intervalroot
#      rts=uroot()
#      r2=rsq(rts$root,data=spruce.df)
#      abline(v=seq(floor(min(d)), ceiling(max(d)), by=1), lwd=0.5,col="pink")
#      abline(v=rts$root,h=rsq(rts$root,data=spruce.df))
#      text(rts$root, r2*1.05,paste("knot is:", rts$root))
#      axis(3, round(rts$root,4))
#      axis(4, round(r2,4),col="Red")

#   })

   output$R2 <- renderPlot({
      dsmoothx = seq(min(d),max(d),length=100)
      dsmoothy = seq(min(d),max(d),length=100)
      points = as.matrix(expand.grid(dsmoothx,dsmoothy))
      r2 = map2_dbl(points[,1],points[,2],~rsqmult(.x,.y,data=spruce.df))
      #r2 = data.frame(r2)
      #r2 %>% filter(r2>.7)
      s3d <- scatterplot3d(points[,1],points[,2],r2, highlight.3d=TRUE,
                           angle=55, scale.y=0.7, pch=16,xlab="xk1",ylab="xk2",
                              main="R^2 vs xk1 vs xk2", zlim=c(.7,.9))

   })

   output$knots<-renderTable({

      argmat = matrix(rep(102,0),nrow=1,ncol=102,byrow=TRUE)
      j=0
      dsmoothx = seq(min(d),max(d),length=100)
      dsmoothy = seq(min(d),max(d),length=100)
      points = as.matrix(expand.grid(dsmoothx,dsmoothy))
      r2 = map2_dbl(points[,1],points[,2],~rsqmult(.x,.y,data=spruce.df))
      modr2 = ceiling(r2*100)/100
      for(i in 1:length(modr2)){
         if(modr2[i] == max(modr2)){
            j=j+1
            argmax = i
            argmat[1,j] = argmax
         }
      }
      crit = which.max(r2[argmat])
      crit_index = argmat[crit]

      xk1 = points[crit_index,1]
      xk2 = points[crit_index,2]
      df <- matrix(c(xk1, xk2, max(r2)),nr=1, nc=3, byrow=TRUE)
      colnames(df) = c("xk1 Critical","xk2 Critical","R2 Max")
      df <- data.frame(df)
      format.data.frame(df,7,7)

   })

   output$tab <- renderTable({

      argmat = matrix(rep(102,0),nrow=1,ncol=102,byrow=TRUE)
      j=0
      dsmoothx = seq(min(d),max(d),length=10)
      dsmoothy = seq(min(d),max(d),length=10)
      points = as.matrix(expand.grid(dsmoothx,dsmoothy))
      r2 = map2_dbl(points[,1],points[,2],~rsqmult(.x,.y,data=spruce.df))
      modr2 = ceiling(r2*100)/100
      for(i in 1:length(modr2)){
         if(modr2[i] == max(modr2)){
            j=j+1
            argmax = i
            argmat[1,j] = argmax
         }
      }
      crit = which.max(r2[argmat])
      crit_index = argmat[crit]

      xk1 = points[crit_index,1]
      xk2 = points[crit_index,2]
      delta = input$delta
      perturb1 = seq(xk1-delta,xk1+delta,length=10)
      perturb2 = seq(xk2-delta,xk2+delta,length=10)
      points_approx = as.matrix(expand.grid(perturb1,perturb2))
      approx_max = as.matrix(map2_dbl(points_approx[,1],points_approx[,2],
                                          ~rsqmult(.x,.y,data=spruce.df)))
      colnames(approx_max)=c("Approximate Maximum R2")
      df <- data.frame(approx_max)
      colnames(points_approx) = c("Approximate xk1 Critical","Approximate xk2 Critical")
      x=format.data.frame(df,7,7)
      cbind(points_approx,x)



})

}

# Run the application
shinyApp(ui = ui, server = server)