在R中基于分数多项式GLM给定概率求解协变量值

我基于分数多项式广义线性模型（GLM），希望找到能输出指定概率的协变量值。

模拟数据的R代码 #

# 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)

模型拟合代码 #

library(mfp)

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

当前近似实现方法 #

我需要针对sex的每个水平（"Male"或"Female"），找到使概率等于指定值frontier <- 0.8的duration值。目前只能通过以下近似方法实现：

pred_duration <- seq(10, 20, by = 0.1)
pred <- data.frame(expand.grid(duration = pred_duration,
                       sex = sex),
           p = predict(model1, 
                       newdata = expand.grid(duration = pred_duration,
                                            sex = sex),
                       type = "response"))

pred[which(pred$p > 0.8), ] %>%
  group_by(sex) %>%
  summarize(min(duration))

但我需要精确解。

内容的提问来源于stack exchange，提问作者jackahall