gamlss-dev / gamlss

gamlss: Generalized Additive Models for Location Scale and Shape
https://CRAN.R-project.org/package=gamlss
10 stars 4 forks source link

predict and shiny #14

Open lemaitreh opened 3 days ago

lemaitreh commented 3 days ago

I'm trying to use the predict function within a shiny app and it seems to have a problem with the newdata option:

RStudio 2021.09.0 Build 351 R 4.4.1 gamlss 5.4.22

Not working (predict.gamlss and newdata):

# Load necessary libraries
library(shiny)
library(gamlss)
library(ggplot2)

# dataset
data(abdom)

# Define UI for the application
ui <- fluidPage(
  titlePanel("GAMLSS Model with Random Effect"),
  sidebarLayout(
    sidebarPanel(
      selectInput("variable", "Choose a variable:", 
                  choices = names(abdom), 
                  selected = "global_mean")
    ),
    mainPanel(
      plotOutput("plot")
    )
  )
)

# Define server logic
server <- function(input, output) {

  output$plot <- renderPlot({
    # Extract the chosen variable
    variable <- input$variable

    # Create a data frame
    holder <- data.frame(y = abdom[[variable]], x = abdom$x)
    holder <- holder[order(holder$x),]

    # Full model
    model <- gamlss(
      y ~ poly(x, 2),
      sigma.formula = ~ poly(x, 2),
      nu.formula = ~ poly(x, 1),
      family = GG,
      data = holder,
      trace = TRUE,
      control = gamlss.control(n.cyc = 100)
    )

    holder$prediction <- predict(model, what = "mu", type = "response", newdata = holder)

    # Plot the data
    ggplot() +
      geom_point(data = holder, aes(x = x, y = y, ), alpha = 0.5) +
      labs(title = "GAMLSS Model",
           x = "X",
           y = "Y") +
      theme_minimal()
  })
}

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

Error:

Avis : Error in eval: objet 'holder' introuvable
  174: eval
  173: eval
  172: predict.gamlss
  170: renderPlot [#22]
  168: func
  128: drawPlot
  114: <reactive:plotObj>
   98: drawReactive
   85: renderFunc
   84: output$plot
    3: runApp
    2: print.shiny.appobj
    1: <Anonymous>

Working (gamlss, no newdata):

# Load necessary libraries
library(shiny)
library(gamlss)
library(ggplot2)

# dataset
data(abdom)

# Define UI for the application
ui <- fluidPage(
  titlePanel("GAMLSS Model with Random Effect"),
  sidebarLayout(
    sidebarPanel(
      selectInput("variable", "Choose a variable:", 
                  choices = names(abdom), 
                  selected = "global_mean")
    ),
    mainPanel(
      plotOutput("plot")
    )
  )
)

# Define server logic
server <- function(input, output) {

  output$plot <- renderPlot({
    # Extract the chosen variable
    variable <- input$variable

    # Create a data frame
    holder <- data.frame(y = abdom[[variable]], x = abdom$x)
    holder <- holder[order(holder$x),]

    # Full model
    model <- gamlss(
      y ~ poly(x, 2),
      sigma.formula = ~ poly(x, 2),
      nu.formula = ~ poly(x, 1),
      family = GG,
      data = holder,
      trace = TRUE,
      control = gamlss.control(n.cyc = 100)
    )

    holder$prediction <- predict(model, what = "mu", type = "response")

    # Plot the data
    ggplot() +
      geom_point(data = holder, aes(x = x, y = y, ), alpha = 0.5) +
      labs(title = "GAMLSS Model",
           x = "X",
           y = "Y") +
      theme_minimal()
  })
}

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

Working (predict.lm, newdata):

# Load necessary libraries
library(shiny)
library(gamlss)
library(ggplot2)

# dataset
data(abdom)

# Define UI for the application
ui <- fluidPage(
  titlePanel("GAMLSS Model with Random Effect"),
  sidebarLayout(
    sidebarPanel(
      selectInput("variable", "Choose a variable:", 
                  choices = names(abdom), 
                  selected = "global_mean")
    ),
    mainPanel(
      plotOutput("plot")
    )
  )
)

# Define server logic
server <- function(input, output) {

  output$plot <- renderPlot({
    # Extract the chosen variable
    variable <- input$variable

    # Create a data frame
    holder <- data.frame(y = abdom[[variable]], x = abdom$x)
    holder <- holder[order(holder$x),]

    # Full model
    model <- lm(
      y ~ poly(x, 2),
      data = holder)

    holder$prediction <- predict(model, what = "mu", type = "response", newdata = holder)

    # Plot the data
    ggplot() +
      geom_point(data = holder, aes(x = x, y = y, ), alpha = 0.5) +
      labs(title = "GAMLSS Model",
           x = "X",
           y = "Y") +
      theme_minimal()
  })
}

# Run the application 
shinyApp(ui = ui, server = server)
lemaitreh commented 3 days ago

temporary fix with the use of global assignment operator:

 holder <<- data.frame(y = abdom[[variable]], x = abdom$x)
zeileis commented 2 days ago

Thanks for the report. I can confirm the bug. The newdata is actually found without any problem. But for some reason the predict() method tries to re-evaluate the Call$data from the fitted model and it does so in the wrong environment. Not sure why this is neeeded in the first place. A minimal reproducible example is:

library("gamlss")
pred <- function() {
  d <- data.frame(y = sin(1:10), x = 1:10)
  m <- gamlss(y ~ x, data = d, family = NO)
  predict(m, newdata = d)
}
pred()
## GAMLSS-RS iteration 1: Global Deviance = 20.7488 
## GAMLSS-RS iteration 2: Global Deviance = 20.7488 
## Error in eval(Call$data) : object 'd' not found

Remark: Using the successor package gamlss2 (not yet on CRAN) the same code works. But the predict() method works a bit differently.