Support the bivariate normal distribution (Solves #12)

DESCRIPTION

@@ -12,7 +12,8 @@ Imports:
         ggplot2 (>= 2.2.1),
         rhandsontable (>= 0.3.4),
         tidyr (>= 0.6.3),
-        rstudioapi (>= 0.6)
+        rstudioapi (>= 0.6),
+        plotly (>= 4.7.0)
 Description: bamlss.vis calls a Shiny app where one can compare different predictions depending on varying covariate values.
 License: GPL-2
 LazyData: TRUE

NAMESPACE

@@ -9,4 +9,11 @@ import(ggplot2)
 import(rhandsontable)
 import(rstudioapi)
 import(shiny)
+importFrom(mvtnorm,pmvnorm)
+importFrom(plotly,add_surface)
+importFrom(plotly,ggplotly)
+importFrom(plotly,plot_ly)
+importFrom(plotly,plotlyOutput)
+importFrom(plotly,plotly_empty)
+importFrom(plotly,renderPlotly)
 importFrom(tidyr,gather)

R/Untitled.R

@@ -0,0 +1,24 @@
+library(shiny)
+library(plotly)
+
+ui <- fluidPage(
+  plotlyOutput("plotly"),
+  plotOutput("plot")
+)
+
+server <- function(input, output) {
+
+  # I need to be able to put NULL in here or anything so that
+  # there is no output and also no error message.
+  output$plotly <- renderPlotly({
+    NULL
+  })
+
+  # This works and sends no error message
+  output$plot <- renderPlot({
+    NULL
+  })
+
+}
+
+shinyApp(ui, server)

R/is.continuous.R → R/is.stuff.R

@@ -12,3 +12,14 @@ is.continuous <- function(family) {
   else
     stop("Family not implemented")
 }
+
+#' Internal: Function to check whether the modeled response is bivariate
+#'
+#'
+
+is.2d <- function(family, links) {
+  if (length(links) == 5 & family == ".mvnorm")
+    return(TRUE)
+  else
+    return(FALSE)
+}

R/plot_dist.R

@@ -48,10 +48,13 @@ plot_dist <- function(model, predictions, palette = "default",
     stop("There is no cdf for the multinomial distribution!")
 
   # Get plot limits
-  lims <- bamlss.vis:::limits(p_m, fam)
+  if (!is.2d(fam, fam_gen$links))
+    lims <- bamlss.vis:::limits(p_m, fam)
 
   # Different plots depending on type of distribution
-  if (bamlss.vis:::is.continuous(fam))
+  if (bamlss.vis:::is.2d(fam, fam_gen$links))
+    plot <- pdfcdf_2d(p_m, model, type)
+  else if (bamlss.vis:::is.continuous(fam))
     plot <- pdfcdf_continuous(lims, funs_list, type, p_m, palette)
   else if (!bamlss.vis:::is.continuous(fam))
     plot <- pdfcdf_discrete(p_m, palette, fam, type, model)
@@ -60,38 +63,6 @@ plot_dist <- function(model, predictions, palette = "default",
   return(plot)
 }
 
-#' Internal: Get plot limits for a predicted distribution.
-#'
-#' Returns a data.frame
-
-limits <- function(predictions, family, times_sd = 3) {
-  if (bamlss.vis:::is.continuous(family)) {
-    if (family == "beta") {  # beta can only be from 0 to 1
-      return(data.frame(x = c(0, 1)))
-    } else if (family == "Generalized Pareto") {
-      return(data.frame(x = c(0, 5)))
-    } else {
-      # Moments
-      moments <- as.data.frame(bamlss.vis:::moments(predictions, family))
-
-      # Limits for each 2 moments
-      lims <- apply(moments, 1, function(x)
-        return(c(x[1] - times_sd * sqrt(x[2]),
-                 x[1] + times_sd * sqrt(x[2]))))
-
-      return(data.frame(x = c(min(lims), max(lims))))
-    }
-  } else if (!bamlss.vis:::is.continuous(family)) {
-    if (family == "binomial") {
-      return(data.frame(x = c(0, 1)))
-    } else if (family == "poisson") {
-      # poisson
-    } else if (family == "multinomial") {
-      # multinomial here
-    }
-  }
-}
-
 #' Internal: Create the pdf/cdf for continuous covariates
 #'
 #' Returns a plot
@@ -193,55 +164,29 @@ pdfcdf_discrete <- function(p_m, palette, family, type, model) {
   return(ground)
 }
 
-#' Internal: Transform discrete predictions into a usable df
+#' Internal: Create 3D pdf/cdf plot
 #'
-#' @importFrom tidyr gather
-
-disc_trans <- function(predictions, family, type, model) {
-  if (family == "binomial") {
-    if (type == "pdf") {
-      predictions$pi_inv <- 1 - predictions$pi
-      predictions$rownames <- row.names(predictions)
-      colnames(predictions) <- c("0", "1", "rownames")
-      tf_df <- gather(predictions, "type", "value", -rownames)
-    } else if (type == "cdf") {
-      predictions$pi_inv <- 1
-      predictions$rownames <- row.names(predictions)
-      colnames(predictions) <- c("0", "1", "rownames")
-      tf_df <- gather(predictions, "type", "value", -rownames)
-      tf_df <- rbind(tf_df, data.frame(rownames = unique(tf_df$rownames),
-                                       type = rep(-1e-100, (nrow(tf_df)/ 2)), # this is because starting point has to be left by just a little margin for plot...
-                                       value = rep(0, (nrow(tf_df)/ 2))))
-      tf_df$type <- as.numeric(tf_df$type)
-    }
-  } else if (family == "poisson") {
-    if (type == "pdf"){
-      limits <- 0:((max(predictions$lambda)*2) + 3) # what lim should preds be?
-      tf_df <- apply(predictions, 1, FUN = function(x) return(dpois(limits, x)))
-      tf_df <- cbind(tf_df, limits)
-      colnames(tf_df) <- c(row.names(predictions), "type")
-      tf_df <- gather(as.data.frame(tf_df), key = "rownames", "value", -type)
-    } else if (type == "cdf") {
-      limits <- 0:((max(predictions$lambda)*2) + 3) # what lim should preds be?
-      tf_df <- apply(predictions, 1, FUN = function(x) return(ppois(limits, x)))
-      tf_df <- cbind(tf_df, limits)
-      colnames(tf_df) <- c(row.names(predictions), "type")
-      tf_df <- gather(as.data.frame(tf_df), key = "rownames", "value", -type)
-      tf_df <- rbind(data.frame(type = 0, rownames = row.names(predictions),
-                     value = -1e-100), # this is because starting point has to be left by just a little margin for plot...
-                     tf_df)
-    }
-  } else if (family == "multinomial") {
-    if (type == "pdf") {
-      tf_df_start <- mult_trans(predictions, model)
-      tf_df <- tf_df_start
-      tf_df$rownames <- row.names(tf_df)
-      tf_df <- gather(tf_df, "type", "value", -rownames)
-      tf_df$type <- factor(tf_df$type, labels = colnames(tf_df_start))
-    }
+#' @importFrom plotly plot_ly add_surface
+
+pdfcdf_2d <- function(p_m, model, type) {
+  if (type == "pdf") {
+    density_f <- mvnorm_bamlss(k = 2)$d
+    p <- p_m[nrow(p_m), ] # only last predicton will be evaluated
+    xval <- seq(p$mu1 - 3 * p$sigma1, p$mu1 + 3 * p$sigma1, length.out = 100)
+    yval <- seq(p$mu2 - 3 * p$sigma2, p$mu2 + 3 * p$sigma2, length.out = 100)
+    z <- matrix(0, ncol = 100, nrow = 100)
+    for (i in 1:100)
+      for (j in 1:100)
+        z[i, j] <- density_f(cbind(xval[i], yval[j]), par = p)
+    return(plot_ly(x = xval, y = yval, z = z) %>% add_surface())
+  } else if (type == "cdf") {
+    p <- p_m[nrow(p_m), ] # only last prediction will be evaluated
+    xval <- seq(p$mu1 - 3 * p$sigma1, p$mu1 + 3 * p$sigma1, length.out = 100)
+    yval <- seq(p$mu2 - 3 * p$sigma2, p$mu2 + 3 * p$sigma2, length.out = 100)
+    z <- matrix(0, ncol = 100, nrow = 100)
+    for (i in 1:100)
+      for (j in 1:100)
+        z[i, j] <- real_pmvnorm(c(xval[i], yval[j]), par = p)
+    return(plot_ly(x = xval, y = yval, z = z) %>% add_surface())
   }
-  return(tf_df)
 }
-
-
-

R/plot_related_helper.R

@@ -0,0 +1,84 @@
+#' Internal: Get plot limits for a predicted distribution.
+#'
+#' Returns a data.frame
+
+limits <- function(predictions, family, times_sd = 3) {
+  if (family == ".mvnorm")
+    return(NULL)
+   else if (bamlss.vis:::is.continuous(family)) {
+    if (family == "beta") {  # beta can only be from 0 to 1
+      return(data.frame(x = c(0, 1)))
+    } else if (family == "Generalized Pareto") {
+      return(data.frame(x = c(0, 5)))
+    } else {
+      # Moments
+      moments <- as.data.frame(bamlss.vis:::moments(predictions, family))
+
+      # Limits for each 2 moments
+      lims <- apply(moments, 1, function(x)
+        return(c(x[1] - times_sd * sqrt(x[2]),
+                 x[1] + times_sd * sqrt(x[2]))))
+
+      return(data.frame(x = c(min(lims), max(lims))))
+    }
+  } else if (!bamlss.vis:::is.continuous(family)) {
+    if (family == "binomial") {
+      return(data.frame(x = c(0, 1)))
+    } else if (family == "poisson") {
+      # poisson
+    } else if (family == "multinomial") {
+      # multinomial here
+    }
+  }
+}
+
+
+#' Internal: Transform discrete predictions into a usable df
+#'
+#' @importFrom tidyr gather
+
+disc_trans <- function(predictions, family, type, model) {
+  if (family == "binomial") {
+    if (type == "pdf") {
+      predictions$pi_inv <- 1 - predictions$pi
+      predictions$rownames <- row.names(predictions)
+      colnames(predictions) <- c("0", "1", "rownames")
+      tf_df <- gather(predictions, "type", "value", -rownames)
+    } else if (type == "cdf") {
+      predictions$pi_inv <- 1
+      predictions$rownames <- row.names(predictions)
+      colnames(predictions) <- c("0", "1", "rownames")
+      tf_df <- gather(predictions, "type", "value", -rownames)
+      tf_df <- rbind(tf_df, data.frame(rownames = unique(tf_df$rownames),
+                                       type = rep(-1e-100, (nrow(tf_df)/ 2)), # this is because starting point has to be left by just a little margin for plot...
+                                       value = rep(0, (nrow(tf_df)/ 2))))
+      tf_df$type <- as.numeric(tf_df$type)
+    }
+  } else if (family == "poisson") {
+    if (type == "pdf"){
+      limits <- 0:((max(predictions$lambda)*2) + 3) # what lim should preds be?
+      tf_df <- apply(predictions, 1, FUN = function(x) return(dpois(limits, x)))
+      tf_df <- cbind(tf_df, limits)
+      colnames(tf_df) <- c(row.names(predictions), "type")
+      tf_df <- gather(as.data.frame(tf_df), key = "rownames", "value", -type)
+    } else if (type == "cdf") {
+      limits <- 0:((max(predictions$lambda)*2) + 3) # what lim should preds be?
+      tf_df <- apply(predictions, 1, FUN = function(x) return(ppois(limits, x)))
+      tf_df <- cbind(tf_df, limits)
+      colnames(tf_df) <- c(row.names(predictions), "type")
+      tf_df <- gather(as.data.frame(tf_df), key = "rownames", "value", -type)
+      tf_df <- rbind(data.frame(type = 0, rownames = row.names(predictions),
+                                value = -1e-100), # this is because starting point has to be left by just a little margin for plot...
+                     tf_df)
+    }
+  } else if (family == "multinomial") {
+    if (type == "pdf") {
+      tf_df_start <- mult_trans(predictions, model)
+      tf_df <- tf_df_start
+      tf_df$rownames <- row.names(tf_df)
+      tf_df <- gather(tf_df, "type", "value", -rownames)
+      tf_df$type <- factor(tf_df$type, labels = colnames(tf_df_start))
+    }
+  }
+  return(tf_df)
+}

R/real_pmvnorm.R

@@ -0,0 +1,9 @@
+#' Internal: Real cumulative density function of multivariate normal dist
+#'
+#' @importFrom mvtnorm pmvnorm
+
+real_pmvnorm <- function(y, par) {
+  return(pmvnorm(upper = y, mean = c(par$mu1, par$mu2),
+                 sigma = matrix(c(par$sigma1^2, 0,
+                                  0, par$sigma2^2), nrow = 2)))
+}

R/vis.R

@@ -5,6 +5,7 @@
 #' @import shiny
 #' @import rhandsontable
 #' @import rstudioapi
+#' @importFrom plotly renderPlotly plotlyOutput ggplotly plotly_empty
 #' @export
 
 ### --- Shiny App --- ###
@@ -35,7 +36,7 @@ vis <- function() {
                         #verbatimTextOutput("testprint"))
                         fluidRow(
                           column(width = 9,
-                                 plotOutput("dist_plot")),
+                                 uiOutput("condition_plot")),
                           column(width = 3, br(),
                                  uiOutput("plotbar"))
                         )
@@ -81,10 +82,24 @@ vis <- function() {
 
     # Reactive model data
     m_data <- reactive({
-      if(!is.null(m()))
+      if (!is.null(m()))
         m()$model.frame
     })
 
+    # Reactive model family
+    fam <- reactive({
+      if (!is.null(m()))
+        family(m())
+    })
+
+    # Got Model and data?
+    gmad <- reactive({
+      if (!is.null(m()) & !is.null(pred$data))
+        TRUE
+      else
+        FALSE
+    })
+
     ## --- Overview tab --- ##
 
     # Equations Output
@@ -248,15 +263,42 @@ vis <- function() {
 
     ## --- Plot Tab --- ##
 
-    ## Plot is rendered here
-    output$dist_plot <- renderPlot({
-      if (!is.null(pred$data))
-        plot_dist(m(), cur_pred(), palette = input$pal_choices,
-                  type = input$type_choices)
+    ## PLotly is rendered here, condition is checked with conditionalPanel
+    output$plotly <- renderPlotly({
+      if (gmad()) {
+        if (is.2d(fam()$family, fam()$links)) {
+          plot_dist(m(), cur_pred(), palette = input$pal_choices,
+                    type = input$type_choices)
+        }
+      } else {
+        plotly_empty()
+      }
+    })
+
+    ## PLotly is rendered here, condition is checked with conditionalPanel
+    output$plot <- renderPlot({
+      if (gmad())
+        if (!is.2d(fam()$family, fam()$links))
+          plot_dist(m(), cur_pred(), palette = input$pal_choices,
+                    type = input$type_choices)
       else
         NULL
     })
 
+    ## The Plot Ui element itself is rendered here
+    ## It checks the conditions for plot and then decides if plotly or plot
+    output$condition_plot <- renderUI({
+      if (gmad()) {
+        if (is.2d(fam()$family, fam()$links)) {
+          plotlyOutput("plotly")
+        } else {
+          plotOutput("plot")
+        }
+      } else {
+        NULL
+      }
+    })
+
     ## Color Choices / pdf/cdf choice are rendered here
     output$plotbar <- renderUI({
       if (!is.null(m()) & any(input$pillpanel == 5, input$pillpanel == 6)) {