Hypothesis test for the difference in Logistic regression probabilities (R glm) 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')
                 )

I get a standard error of the probability output for one observation:
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))

So, I am able to get probability p1, and standard error SE1 respective to a feature vector X1 of observation A (a specific person doing something at a specific moment in time). Other observations (Xi) used to fit the logistic model are randomly sampled and independent.
Then, suppose I have a new feature vector X2 with respect to observation A, which has a physical meaning of how the person behaved at another moment of time, which does not overlap with the first time period. I assume here that X1 and X2 are not necessarily dependent, so they can be arbitrarily different.
I run the new vector of features through my trained model and get probability p2 and standard error of prediction SE2. 
I wanto to test if p1 is significatly different from p2. Z = (p1 - p2) / pooled_se.
Question: A confusion arises when 1) I decide how to get a pooled standard deviation, and 1.a) what are the degrees of freedom? Are these the number of observations used in model fitting, or are these equal to 1, since I run a hypothesis test for just one observation?
 A: I think I got an answer for equal variances and equal sample size:
xx <- rnorm(100, 0, 2) #Population SD = 2
yy <- rnorm(100, 0, 2) #Population SD = 2

# xx <- rnorm(100, 0, 1) #Population SD = 1
# yy <- rnorm(100, 0, 2) #Population SD = 2



SEp <- 2 * sqrt(1/100 + 1/100) # pooled SE with equal population variances

SEp <- sqrt(1^2/100 + 2^2/100) #for unknown variances, using sample SDs

# sample-wise standard errors

SEx <- 1/sqrt(100)
SEy <- 1/sqrt(100)


# given SEs only

SEp_alt <- sqrt(SEx^2 + SEy^2)

print(paste0('full calculation = ', SEp, '; ', 'SE only calculation = ', SEp_alt))

So the answer would be; pooled SE equals
sqrt(SEx^2 + SEy^2)

which matches the pooled SE result with both known and unknown population standard deviation, and equal sample sizes.
A: I tried to find the correlations between probabilities returned by logistic regression model. I mean Pearson correlations for two random variables of probabilities. I build 5 000 sets of random feature vectors, assuming that each observation in any set is pairwise connected to the observation in any other set (same entities are under testing). I assume that the model quality is good, and there is a linear separability of classes corrupted my controlled level of Guassian noise.
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'))

I did not find significant correlations.

So I am asking again why you say the p values stemming from one model with different, independent features as an input would be correlated. Are there specific assumptions that should be met to satisfy this?
