Проверка гипотезы о различии в вероятностях из модели логистической регрессии

dummy_data <- data.frame(
     x1 = rnorm(100, 0, 1),
     x2 = rnorm(100, 0, 1),
     y = round(runif(100, 0, 1), 0)
)
 
dummy_glm <- glm(y ~ x1 + x2,
                 data = dummy_data,
                 family = binomial(link = 'logit')
                 )
 
#get var-covar matrix of coefficient estimates
 
X <- as.matrix(cbind(rep(1, times = 100), dummy_data[, c('x1', 'x2')]), ncol = 3)
 
### FOR LOGISTIC RGRESSION
 
pi <- dummy_glm$fitted.values
 
w <- pi * (1 - pi)
 
v <- diag(w, length(w), length(w))
 
XtX_inv <- solve(t(X) %*% v %*% X)
 
# standard error for linear combination (first observation in data frame)
 
C <- c(1, dummy_data[1, 'x1'], dummy_data[1, 'x2'])
 
std_err_manually <- as.numeric(sqrt(t(C) %*% XtX_inv %*% C))

rm(list = ls()); gc()
 
library(data.table)
 
#dummy_data <- read.csv('dummy_data.csv')
 
dummy_data <- data.table(
     x1 = rnorm(100, 0, 1),
     x2 = rnorm(100, 0, 1)
)
dummy_data[, y := ifelse(x1 + x2 + rnorm(1, 0, 0.5) >= 0, 1, 0)]
 
 
dummy_glm <- glm(y ~ x1 + x2 -1,
                 data = dummy_data,
                 family = binomial(link = 'logit')
                 )
 
summary(dummy_glm)
 
# run simulation
iterations <- 5000
fit_list <- list()
 
for (i in 1:iterations){
 
     fit_list[[i]] <- predict(
          object = dummy_glm,
          newdata = dummy_data,
          se.fit = T,
          type = 'response'
     )$fit
 
     # make new independent inputs vars
 
     dummy_data <- data.table(
          x1 = rnorm(100, 0, 1),
          x2 = rnorm(100, 0, 1)
     )
     dummy_data[, y := ifelse(x1 + x2 + rnorm(1, 0, 0.5) >= 0, 1, 0)]
 
}
 
fit_probabilities <- t(do.call(rbind, fit_list))
 
fit_prob_corr_matrix <- cor(fit_probabilities)
 
x <- fit_prob_corr_matrix[lower.tri(fit_prob_corr_matrix, diag = FALSE)]
 
hist(x, breaks = 'fd', main = paste0('density over Pearson Corr Coef after ', iterations, ' iterations'))