EDIT: The tl;dr here would broadly be: given that both frequentist standard errors and a quadratic approximation of a Bayesian joint posterior can be obtained from the square root of the diagonal elements of the inverse of the negative Hessian of the log-likelihood function at the MLE (ie, the FIM), when can I use sumstats from a frequentist model fit for posterior predictive simulation? Sorry for the length, new to this site!
ORIGINAL TEXT:
So I have a question about linear models, specifically in the context of their use in prediction while accommodating inferential uncertainty. Normally, I would just fit the model in question numerically in a Bayesian framework, using samples from the joint posterior to approximate a posterior predictive distribution of whatever output variable given whatever novel input variable.
In this context, though, I don't have access to the original data from the fitted linear model. Instead, I have access to a bunch of p-values & signs of effect, & can approximate the variance of my input variables (to further clarify, this is specifically motivated by the idea of a polygenic score, though I'm equally as interested in the generality of this sorta thing).
Some thought later and I don’t think we're working in the world of Bayesian sufficiency (if that's the right term?) here. However, I do think these values are enough for a quadratic approximation of the log-likelihood surface, so long as there’s not too much non-identifiability and the whole thing’s roughly parabolic, which I think should be asymptotically true for flat models like this. And thus also enough for a quadratic approximation of the posterior distribution. Backing up a step to play around in R, let’s specify and fit a linear model:
## specify generative model parameters
n = 20
a = 3
b = 0.1
## simulate data
x_prob <- 0.23
x = rbinom(n = n, prob = x_prob, size = 2)
y = a + b * x + rnorm(n)
d <- data.frame(x = x, y = y)
d$ys <- (d$y - mean(d$y)) / sd(d$y)
d$xs <- (d$x - mean(d$x)) / sd(d$x)
fit <- lm(y~x, d)
We can go forwards to find the p-value of the coefficient on x just from knowing the Pearson correlation between x and y and the sample size (unsurprising, since Pearson correlations are just standardized regression coefficients):
## from correlation to pval
summary(fit)$coefficients[2,4]
(1 - pnorm(abs(atanh(cor(d$x,d$y))) / (1 / sqrt(n-3))))*2 #using Fisher transf. / approx.
tail_prop <- function(p) min(p, 1-p) * 2
tail_prop(pt(q = cor(d$x,d$y) * sqrt(n-2) / sqrt(1-cor(d$x,d$y)^2), df = n - 2))
cor.test(x, y)$p.value #just to confirm
Inverting this procedure, we can find the correlation between variables from the p-value and a direction:
## from pval to correlation
pval <- summary(fit)$coefficients[2,4]
tanh(qnorm(1-pval/2) / sqrt(n-3)) * sign(summary(lm(y~x, d))$coefficients[2,1]) #using Fisher transf. / approx.
inv_tail_prop_right <- function(p) 1 - p / 2
a = qt(p = inv_tail_prop_right(pval), df = n - 2) / sqrt(n-2)
abs(a / sqrt(a^2 +1)) * sign(summary(fit)$coefficients["x", "Estimate"])
cor(d$x,d$y)
Since we’re interested in the count-specific effect, we can un-standardize the correlation for x, leaving the outcome in units sd(y). Since sd(x) tends to not be reported, we can assume the population is at HW equilibrium and approximate it:
## from pval to partly standardized coefficient
pval <- summary(fit)$coefficients[2,4]
tanh(qnorm(1-pval/2) / sqrt(n-3)) / sqrt(x_prob * (1-x_prob) * 2) * sign(summary(lm(ys~x, d))$coefficients[2,1]) #using Fisher transf. / approx.
summary(lm(ys~x, d))$coefficients[2,1]
a = qt(p = inv_tail_prop_right(pval), df = n - 2) / sqrt(n-2)
abs(a / sqrt(a^2 +1)) * sign(summary(fit)$coefficients["x", "Estimate"]) / sqrt(x_prob * (1-x_prob) * 2)
#last bit is the variance of a binomial random variable, which would assume pop at HW-equilibrium. For exact equality to lm(), divide by sd(x), which we don't know
abs(a / sqrt(a^2 +1)) * sign(summary(fit)$coefficients["x", "Estimate"]) / sd(x)
The standard errors of the MLE do tend to be reported, and I think they can be used to perform a quadratic approximation of the underlying likelihood surface (or joint posterior distribution, if it’s approx. normal). Specifically, I think this is because the Hessian of the likelihood (or more easily, the log-likelihood) can be used to both find standard errors, as well as describe a quadratic approximation of the joint posterior from a flat prior (or any prior, so long as you re-weight & re-scale things):
## can use standard errors for quadratic approx of marginal posterior
library(calculus)
#specify the likelihood
likelihood_string <- sapply(1:n, function(i) paste0("(1 / sig / sqrt(2 * pi) * exp(-1/2 * ((",
d$ys[i], "-(a+b*", d$x[i], "))/sig)^2))"))
likelihood <- parse(text = paste0(likelihood_string, collapse = " * "))
vars <- c(sig = summary(fit)$sigma,
a = summary(fit)$coefficients["(Intercept)","Estimate"],
b = summary(fit)$coefficients["x","Estimate"])
#compute & compare the hessian
hess_mat <- hessian(f = likelihood, var = vars)
summary(fit)$coefficients[,"Std. Error"]
sqrt(diag(solve(-hess_mat / prod(dnorm(x = d$ys,
mean = summary(fit)$coefficients[1,1] + summary(fit)$coefficients[2,1] * d$x,
sd = summary(fit)$sigma)))))[-1]
#now try the log-likelihood
loglikelihood_string <- sapply(1:n, function(i) paste0("log(1 / sig / sqrt(2 * pi) * exp(-1/2 * ((",
d$ys[i], "-(a+b*", d$x[i], "))/sig)^2))"))
loglikelihood <- parse(text = paste0(loglikelihood_string, collapse = " + "))
vars <- c(sig = summary(fit)$sigma,
a = summary(fit)$coefficients["(Intercept)","Estimate"],
b = summary(fit)$coefficients["x","Estimate"])
hess_mat <- hessian(f = loglikelihood, var = vars)
summary(fit)$coefficients[,"Std. Error"]
sqrt(diag(solve(-hess_mat)))[-1]
We can manually confirm that the marginal posterior distribution coming from a flat prior has first and second moments (ie, completely describing a corresponding normal) equal to the OLS estimate and squared standard error:
## manually find first and second moment of marginal posterior
moments_from_grid <- function(x, d, m = c(1:2), ll = F){
#prevent underflow -- gets normalized out later
if(ll){
d <- d - max(d) + 10
d <- exp(d)
}
m1 <- sum(x*d) / sum(d)
m_eval <- setdiff(m, 1)
ms <- sapply(m_eval, function(m_i) sum((x-m1)^m_i*d) / sum(d))
if(1 %in% m){return(unlist(c(m1, ms)))}else{return(ms)}
}
b_range <- -500:500/100
likelihood_surface_b <- data.frame(b = b_range, ll = sapply(b_range, function(bi) sum(dnorm(x = d$ys, mean = 0 + d$xs*bi, sd = summary(lm(ys~xs, d))$sigma, log = T))))
# likelihood_surface_b$ll <- likelihood_surface_b$ll - max(likelihood_surface_b$ll) + 3
plot(likelihood_surface_b$b, exp(likelihood_surface_b$ll), type = "l")
c(moments_from_grid(likelihood_surface_b$b, likelihood_surface_b$ll, 1, T),
sqrt(moments_from_grid(likelihood_surface_b$b, likelihood_surface_b$ll, 2, T)))
summary(lm(ys~xs, d))$coefficients[2,1:2]
If the model fit also reports the coefficient of determination, I think we can use it to approximate the uncertainty left in the outcome conditional on the predictors. If they don’t, we can still estimate it from the reported p-values, similar to how we would estimate SNP-based heritability:
(1 - summary(lm(ys~0+xs, d))$sigma)*2
summary(lm(ys~0+xs, d))$r.squared
So to sum, does this all seem reasonable? Can I (making a few other simplifying assumptions, e.g. independence between inputs) compute posterior predictive distributions in this manner?
(I think one of my main concerns pertain to the case where the underlying model’s non-identifiable, and the true Hessian / FIM has tons of non-trivial off-diagonal elements and is uninvertible etc. Then, using standard errors for the quadratic approximation would be invalid, I think)
