3
$\begingroup$

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)

$\endgroup$
5
  • 4
    $\begingroup$ This post is so long that I went straight to the end in my search for a clear question. "Does all this seem reasonable" looks like a very broad unfocused question and I just don't want to have to read through pages of preliminary material to determine what the question(s) might be. I suspect most readers will be in a similar position. If you can find a way to abbreviate your post, you are more likely to get a useful answer. $\endgroup$
    – whuber
    Aug 21 at 14:52
  • $\begingroup$ Gotcha -- sorry for the lengthy preamble, I'd wanted to describe my broader thought process in framing the question, as well as validate my intuitions with some basic numerical results. The tl;dr 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? Will edit this in above! $\endgroup$ Aug 21 at 15:55
  • 4
    $\begingroup$ That's a good idea. I believe many of us appreciate the attention to detail you provide. My suggestion is to preface your question with a summary much like that in your comment and follow that with details for those who wish to read further. $\endgroup$
    – whuber
    Aug 21 at 17:57
  • $\begingroup$ Could you also edit the "The tl;dr here" sentence? $\endgroup$
    – Xi'an
    Aug 25 at 11:35
  • $\begingroup$ Sure -- how so? $\endgroup$ Aug 25 at 14:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.