I have a quadratic term in a GLM and I am interested in the vertex value (+ the standard error and confidence interval of the vertex) of the quadratic term. To my knowledge, there is no automatic function for this purpose in R, and I am unable to calculate this manually, as I do not have sufficient knowledge of statistics and R. Is there anybody who could help and build the code to derive the values of interest?
GLM <- glm(Y ~ X1 + I(X1^2) + X2 + X3 + X4 + X5, family="binomial", data=DF)
Warnings aside since I am not used to simulating data for these types of problems. Here is a vertex with CI based on a Bonferroni CI of the parameters.
Since we need two parameters, we need to adjust the product of the two CIs such that x^2=1-alpha -> x=(1-alpha)^(1/2)
set.seed(1)
library(dplyr, warn.conflicts = F)
library(ggplot2, warn.conflicts = F)
logit_trans <- function(x) 1/(1/exp(x)+1)
sim_data <- tibble(
x = rep(seq(0,20,.1), 10),
y = round(logit_trans(20 - 8*x + .4*x^2 + rnorm(2010, 0, 1)))
)
model <- glm(
y ~ x + I(x^2),
binomial(),
sim_data
)
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
find_vertex <- function(model_pars){
a <- model_pars[2]
b <- model_pars[3]
if(2*b>0){ #Check second deriv.
-a/(2*b)
} else{ # concave down
NA
}
}
ci_vertex <- function(model){
mean_vertex <- model %>%
coef() %>%
find_vertex()
ci <- model %>%
confint(parm = 2:3, level = (1-.05)^(1/2))
lower_ci <- ci %>%
magrittr::extract(,1) %>%
c(0, .) %>%
find_vertex()
upper_ci <- ci %>%
magrittr::extract(,2) %>%
c(0, .) %>%
find_vertex()
list(
lower_ci = lower_ci,
mean_vertex = mean_vertex,
upper_ci = upper_ci
)
}
ci <- ci_vertex(model)
#> Waiting for profiling to be done...
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
sim_data %>%
mutate(
x2 = x^2,
y_expected = logit_trans(c(matrix(c(rep(1, n()), x, x2), ncol = 3) %*% matrix(coef(model))))
) %>%
ggplot(aes(x, y_expected)) +
geom_point() +
geom_line() +
geom_vline(xintercept = ci$lower_ci, color = 'red', linetype = 'dashed') +
geom_vline(xintercept = ci$upper_ci, color = 'red', linetype = 'dashed') +
geom_vline(xintercept = ci$mean_vertex, color = 'red')
Created on 2022-08-19 by the reprex package (v2.0.1)
The mean vertex is almost exactly 10 as we would expect from the input function for y. For some reason, the upper and lower CI are inverted and the lower one is slightly larger than expected. I am guessing this is due to model unidentifiability.
