simple_regression/server.R at main · SurgicalInformatics/simple_regression · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
library(shiny)

shinyServer(function(input, output) {

  # --------------------------------------------------------------------------
  # Get a set of random data with a fixed true model
  draw.sample <- reactive({
    # This gets called whenever the app is reloaded

    # Hardcode the true relationship
    n.obs = 50
    x <- rnorm(n.obs, 3, 2)
    y <- 2 + x + rnorm(n.obs, 0, 1)

    model.summary <- summary(lm(y ~ x))

    return(list(x=x, y=y, model.summary=model.summary))

  })

  # --------------------------------------------------------------------------
  # Calculate the current values of the model given the inputs
  regression <- reactive({

    # Get shorthand access to the attributes we care about
    data.vals <- draw.sample()
    x <- data.vals$x
    y <- data.vals$y
    a <- input$intercept
    b <- input$slope

    # Give a visual cue when we have the right regression
    if (a == 2 & b == 1) resid.color <- "seagreen" else resid.color <- "firebrick"

    # Calculate the current residuals
    yhat <- input$intercept + x * input$slope
    resid <- y - yhat

    # Calculate the current and optimal residual sum squares
    ss.res <- sum(resid ** 2)
    resid.best <- y - (2 + x)
    ss.res.best <- sum(resid.best ** 2)

    # Compute R^2
    r2 <- 1 - (ss.res / sum((y - mean(y)) ** 2))

    return(list(x=x, y=y, yhat=yhat, a=a, b=b, r2=r2,
                resid=resid, ss.res=ss.res, ss.res.best=ss.res.best,
                resid.color=resid.color))

  })

  #---------------------------------------------------------------------------
  # Plot a scatter of the data and the current model with residuals
  output$reg.plot <- renderPlot({

    # Get the current regression data
    reg.data <- regression()
    a <- reg.data$a
    b <- reg.data$b
    x <- reg.data$x
    y <- reg.data$y
    r2 <- reg.data$r2
    resid <- reg.data$resid

    # Mask data outside the viewport
    mask <- x > 0 & x < 4.5 & y > 0 & y < 8
    x <- x[mask]
    y <- y[mask]
    resid <- resid[mask]


    # Plot the regression line
    plot(c(-4.5, 4.5), c(a + b * -4.5,  a + b * 4.5), type="l", lwd=2,
         bty="n", xlim=c(0, 5), ylim=c(0, 8), xlab="Cups of coffee", ylab="Systolic blood pressure minus 120",
         main="Linear model of systolic blood pressure by coffee consumption")

    # Plot each residual distance
    for (i in 1:length(resid)){
      lines(c(x[i], x[i]), c(y[i], y[i] - resid[i]),
            col=reg.data$resid.color, lwd=1.5)
    }

    # Plot the observations
    points(x, y,  pch=16, col="#444444")

    # Plot the current equation as a legend
    legend(-5, 8, sprintf("y = %.3g + %.3g * x", a, b), lty=1, lwd=2, bty="n")

  })

  #---------------------------------------------------------------------------
  # Plot the current sum squares along with the minumum possible
  output$ss.plot <- renderPlot({

    # Get the current regression data
    reg.data <- regression()
    ss.res <- reg.data$ss.res
    ss.res.best <- reg.data$ss.res.best
    resid.color <- reg.data$resid.color

    # Plot the two points
    plot(ss.res, 1, col=resid.color, cex=2,
         yaxt="n", bty="n", xlim=c(0, 1000),
         ylab="", xlab="", main="Sum of Squares of Residuals")
    points(ss.res.best, 1, pch=4, cex=2)

  })

  #----------------------------------------------------------------------------
  # Plot the current distribution of residuals and the theoretical distribution
  output$resid.plot <- renderPlot({

    # Get the current regression data
    reg.data <- regression()
    resid <- reg.data$resid
    resid <- resid[resid > -5 & resid < 5]

    # Plot a histogram of the residuals
    hist(resid, seq(-5, 5, .5), prob=TRUE, col="#bbbbbb",
         xlim=c(-5, 5), ylim=c(0, dnorm(0) * 1.5),
         yaxt="n", bty="n", ylab="", xlab="", main="Distribution of Residuals")
    rug(resid, lwd=2)

    # Plot a normal density (the expected residual distribtuion)
    curve(dnorm, col=reg.data$resid.color, lwd=2, add=TRUE)

  })

  #---------------------------------------------------------------------------
  # Print the glm() summary of the true model
  output$summary <- renderPrint({

    if (input$summary){
      return(draw.sample()$model.summary)
    }

  })
})