How can I find the maximum output of a function-CodePudding

If I have a GLM, is there any way I can efficiently find the maximum output by changing one covariate and holding the others?

Using my simulated data:

# FUNCTIONS ====================================================================

logit <- function(p){
  x = log(p/(1-p))
  x
}

sigmoid <- function(x){
  p = 1/(1   exp(-x))
  p
}

beta_duration <- function(D, select){
  logit(
    switch(select,
           0.05   0.9 / (1   exp(-2*D   25)),
           0.9 * exp(-exp(-0.5 * (D - 11))),
           0.9 * exp(-exp(-(D - 11))),
           0.9 * exp(-2 * exp(-(D - 9))),
           sigmoid(0.847   0.210 * (D - 10)),
           0.7   0.0015 * (D - 10) ^ 2,
           0.7 - 0.0015 * (D - 10) ^ 2   0.03 * (D - 10)
    )
  )
}

beta_sex <- function(sex, OR = 1){
  ifelse(sex == "Female", -0.5 * log(OR), 0.5 * log(OR))
}

plot_beta_duration <- function(select){
  x <- seq(10, 20, by = 0.01)
  y <- beta_duration(x, select)
  data.frame(x = x,
             y = y) %>%
    ggplot(aes(x = x, y = y))  
    geom_line()  
    ylim(0, 1)
}


# DATA SIMULATION ==============================================================

duration <- c(10, 12, 14, 18, 20)
sex <- factor(c("Female", "Male"))

eta <- function(duration, sex, duration_select, sex_OR, noise_sd){
  beta_sex(sex, sex_OR)   beta_duration(duration, duration_select)   rnorm(length(duration), 0, noise_sd)
}

sim_data <- function(durations_type, sex_OR, noise_sd, p_female, n, seed){
  set.seed(seed)
  data.frame(
    duration = sample(duration, n, TRUE),
    sex = sample(sex, n, TRUE, c(p_female, 1 - p_female))
  ) %>%
    rowwise() %>%
    mutate(eta = eta(duration, sex, durations_type, sex_OR, noise_sd),
           p = sigmoid(eta),
           cured = sample(0:1, 1, prob = c(1 - p, p)))
}

# DATA SIM PARAMETERS
durations_type <- 4 # See beta_duration for functions
sex_OR <- 3 # Odds of cure for male vs female (ref)
noise_sd <- 1
p_female <- 0.7 # proportion of females in the sample
n <- 500 

data <- sim_data(durations_type = 1, # See beta_duration for functions
                 sex_OR = 3, # Odds of cure for male vs female (ref)
                 noise_sd = 1,
                 p_female = 0.7, # proportion of females in the sample
                 n = 500,
                 seed = 21874564)

I am fitting a fractional polynomial GLM:

library(mfp)

model1 <- mfp(cured ~ fp(duration)   sex,
              family = binomial(link = "logit"),
              data = data)
summary(model1)

Given that I am holding sex as constant, is there any way to find the value of duration within a certain range that gives me the highest predicted value? Something less inefficient than:

range <- seq(10, 20, by = 1e-4)
range[which.max(predict(model, type = "response", newdata = data.frame(duration = range, sex = "Male")))]

CodePudding user response：

You can use optimize here. Just create a function which returns a prediction based on the value of duration:

f <- function(x) predict(model1, list(sex = 'Male', duration = x))

And we can find the value of duration which produces the maximum log odds within the range 0-20 by doing:

optimise(f, c(0, 20), maximum = TRUE)$maximum
#> [1] 17.95679