shiny_regression.R

shinyApp(
  options = list(width = '25%', height = '25%'),
  ui = shinyUI(fluidPage(
    sidebarLayout(
      sidebarPanel(
        sliderInput('e', label = 'e', min = 0, max = 100, value = 10, step = 1),
        sliderInput('N', label = 'N', min = 10, max = 500, value = 100, step = 1),
        numericInput('seed', label = 'random seed', value = 1774),
        numericInput('lambda', label = 'penalty (lambda)', value = 0),
        textInput('true', label = 'true function f(x)', value = 'x^2'),
        textInput('fitted', label = 'fitted function g(x)', value = 'x'),
        tableOutput('coefficients'),
        hr(),
        p('Sum of squared errors: ', textOutput('SS'))
        ),
      mainPanel(
        plotOutput('update_plot', height = '400px')
        )
      )
    )
  ),
  
  server = function(input, output) {
    
    # compute the ridge regression weights
    ridge <- function(y, X, lambda) {
      cat(lambda, '\n')
      reg <- diag(ncol(X))
      diag(reg) <- lambda
      solve(t(X) %*% X + reg) %*% t(X) %*% y
    }
      
    # parse the true / fitted function string
    parse_string <- function(expr, x) {
      parsed <- NULL
      formsplit <- strsplit(expr, ' ')[[1]]
      len <- length(formsplit)
      parseform <- function(form) eval(parse(text = form), list('x' = x))
      
      while (is.null(parsed)) {
        parsed <- tryCatch({
          parseform(paste(formsplit, collapse = ' ')) 
        }, error = function(e) {
          len <<- len - 1
          formsplit <<- formsplit[1:len]
          NULL
        })
        
        if (len == 0) {
          return(list('parsed' = 'x', 'formsplit' = 'x'))
        }
      }
        
      list('parsed' = parsed, 'formsplit' = formsplit)
    }
    
    # create the design matrix given the parsed function string
    create_mat <- function(fitted_split, x) {
      prepridge <- fitted_split[fitted_split != '+'] # filter all +
      X <- sapply(1:length(prepridge), function(i) {
        eval(parse(text = prepridge[i], list('x' = x)))
      })
      cbind(1, X)
    }
    
    # run the whole regression based on user input
    run_regression <- function(e, N, true, fitted, seed, lambda) {
      if (true == '') true <- 'x'
      if (fitted == '') fitted <- 'x'
      if (is.na(lambda)) lambda <- 0
      if (is.na(seed)) seed <- 1774
      
      set.seed(seed)
      x <- rnorm(N, mean = 0, sd = 5) # predictor
      error <- rnorm(N, mean = 0, sd = e) # error
      
      true_fun <- parse_string(true, x)
      fitted_fun <- parse_string(fitted, x)
      
      y <- true_fun$parsed + error
      X <- create_mat(fitted_fun$formsplit, x)
      betas <- ridge(y, X, lambda)
      pred <- X %*% betas
      
      list('x' = x, 'y' = y, 'X' = X, 'betas' = betas,
           'pred' = pred, 'true_fun' = true_fun)
    }
    
    update_plot <- function(e, N, true, fitted, seed, lambda) {
      reg <- run_regression(e, N, true, fitted, seed, lambda)
      
      # plotting options
      par(cex.main = 1.2, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0),
        cex.lab = 1.2 , font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
      
      plot(reg$x, reg$y, pch = 21, cex = 1.1, bg = 'grey',
           xlab = '', ylab = '', xaxt = 'n', yaxt = 'n')
      par(new = TRUE)
      
      plot(reg$x, reg$true_fun$parsed, col = 'green', lwd = 2, lty = 2, pch = 20,
           cex = .8, yaxt = 'n', xaxt = 'n', ylab = '', xlab = '')
      
      par(new = TRUE)
      plot(reg$x, reg$pred, col = 'red', lwd = 2, lty = 2, pch = 20,
           main = 'regression', cex = .8, ylab = 'y', xlab = 'x')
    }
    
    output$SS <- renderText({
      reg <- run_regression(input$e, input$N, input$true, input$fitted,
                            input$seed, input$lambda)
      round(sum(abs(reg$y - reg$pred)^2), 3)
    })
    
    output$coefficients <- renderTable({
      reg <- run_regression(input$e, input$N, input$true, input$fitted,
                            input$seed, input$lambda)
      n <- length(reg$betas) - 1
      betas <- data.frame('coef' = paste0('b', 0:n), 'value' = reg$betas)
      betas
    })
    
    output$update_plot <- renderPlot({
      update_plot(input$e, input$N, input$true, input$fitted,
                  input$seed, input$lambda)
    })
  }
)