R/lidar_penalization.R

In fairy1991/coolpackage: Interactive graphics for linear regression

#' Penalized model for LIDAR dataset
#'
#' This function allows you to visualize the effect of different values of lambda in the penalized splines, using the LIDAR dataset from SemiPar package.
#'
#' @keywords penalization
#'
#' @import shiny
#' @import ggplot2
#' @export
#'
#' @examples
#' \dontrun{
#' penalized_splines()
#' }
#'
penalized_splines = function(){

	## this is a temporary data
  samplesize = sample(20:50, 1)
  intercept = runif(1, 0, 1)
  slope = runif(1, -intercept, 1 - intercept)

  data.example = data.frame(x = runif(samplesize, 0, 1)) %>%
    mutate(y = intercept + slope * x + rnorm(samplesize, mean = 0, sd = 0.05)) %>%
    arrange(x) %>%
    mutate(., fitted = fitted(lm(y ~ x, data = .)))

	## definitions for penalized smoothing
  range = data.example$x
  y = data.example$y

  knots <- quantile(range, seq(0, 1, length = 45))[-c(1, 45)]
  X.full = cbind(1, range, sapply(knots, function(k) ((range - k > 0) * (range - k))))

  P = diag(1, dim(X.full)[2], dim(X.full)[2])
  P[1,1] = 0

  ## UI definition
  ui <- fluidPage(
    headerPanel('Penalized model for LIDAR dataset'),
  	sidebarPanel(
    	sliderInput(inputId = "lambdaslider",
      	          label = "Choose a value for lambda",
        	        value = 1, min = -2, max = 7, pre = "10^")
  	),
  	mainPanel(
    	plotOutput("lidarplot")
  	)
  )

  ## Server definition
	server <- function(input, output) {
  	lidar.pred <- reactive({
    	data.frame(range = range,
      	         fitted = X.full %*% (solve(t(X.full) %*% X.full + {10^input$lambdaslider} * P) %*% t(X.full) %*% y))
  	})

		output$lidarplot <- renderPlot({
			p <- ggplot(lidar, aes(x = range, y = logratio)) + geom_point() + geom_path(data = lidar.pred(), aes(x = range, y = fitted), color = "blue")
			print(p)
		})
	}

	## run app
	shinyApp(ui = ui, server = server)

}