Survey weighted regression

Very helpful insights, @Guido_Biele, @maxbiostat, @bgoodri.

How could " the …inclusion probabilities [be used] to exponentiate the likelihood contribution of each sampled unit" in a stan program?

Thanks in advance.

In log space that would be

 for(i in 1:N){
    target +=  normal_lpdf(y[i] | yHat[i], sigma) * weights[i];
  }

as pointed out above.

1 Like

Thanks, @maxbiostat. That is very helpful.

It isn’t Bayesian. The W_i people in the population that have the same background characteristics as the i-th person in the sample do not all have y_i as their outcome but multiplying the log-likelihood for the i-th person by W_i assumes so. I would rather post-stratify predictions from a real posterior distribution than work with something that isn’t a real posterior distribution.

2 Likes

Thanks for clarifying! My confusion about how weights are treated by rstanarm stems from the fact that ?stan_lm says that they’re treated the “same as lm”.

I don’t think lm can be taking an approach analogous to normal_lpdf(y[i] | yHat[i], sigma) * weights[i] though, because lm weights are invariant to scale:

> y <- rnorm(100)
> w <- runif(100)
> fit_w   <- lm(y ~ 1, weights = w)
> fit_10w <- lm(y ~ 1, weights = 10*w)
> SE <- function(fit) summary(fit)$coef['(Intercept)', 'Std. Error'] 
> SE(fit_w) == SE(fit_10w)
[1] TRUE

Perhaps it would be clearest to say that weights in rstanarm are replication counts that treat each observation as one or more real observations (analogous to Stata’s fweights).

This would help clarify that scale very much does matter for rstanarm weights – the sum of the weights is equal to the number of observations in the original dataset.

1 Like

The stan_glm and stan_lm functions do things a little differently. For most models, they are interpreted as frequency weights. For stan_lm, they can be like the weights for generalized least squares.

@Guido_Biele. Thanks for this very helpful approach. I have two questions:

(1) How could we make the relationship between the weight and the variance explicit in the model below?

The reason is that we have two types of inverse variance weights:

weight_RandomEffects = 1/(tau^2 + sigma^2)
weight_FixedEffects = 1/sigma^2

(2) Could it be correct to implement directly in a stan program the relevant equations from here, here, here, or elsewhere by adding a transformed parameters block to the model below?

For example,

transformed parameters {
...
vector[N] tau;  
tau = sqrt(sum(w^2[i][(y[i] - y_mu)^2 - sigma^2[i]])/sum(w^2[i]))
...
}

Thanks in advance.

I can look a bit more into this on Monday.
Generally, you can use weights in Stan as you can do it in when computing maximum likelihood estimates. (Though not everyone agrees one should)
The reason I’m am hesitant in my answer is, that I am unsure about what the sigma is in the inverse variance weights. I assume it’s the sigma of the effect sizes, and I mm addition one estimates an error variance, but I can’t tell without access to the paper. (I am also not sure what tau is).

As an aside, you can use brms to estimate meta analysis models. See for example here: https://vuorre.netlify.com/post/2017/01/19/better-forest-plots-from-meta-analytic-models-estimated-with-brms/

Great, @Guido_Biele. Here, we have two sources of variability in the effect sizes of the primary studies:

tau^2 = between-studies variance
sigma^2 = within-study variance

tau can be calculated using the formula shown above in transformed parameters block.

sigma is calculated as follows:

These variances are then used to estimate the inverse variance weights:

Thanks in advance.

Hi Bob,

In a 2017 post, you write

The closest the manual comes is a section on “Exploiting sufficient statistics”. I haven’t added anyting on weighted regression since our regression experts, Ben Goodrich and Andrew Gelman, don’t like the weightings (other than those based on sufficient stats) because the resulting model isn’t properly Bayesian in that there’s no generative process for the weights.

This makes very good sense to me and was just wondering if you could point to a specific paper where Goodrich and Gelman spell out these issues.

Thanks in advance,

David

That would be Gelman (2007)


see also the discussion, rejoinder, and citations for that issue.

1 Like

Thanks!

David

Hi,
I have fitted following logistic regression model with survey weights. I was able to run this model without any error. But it was very slow and it had all sort of warnings including divergent transitions, large R-hat values etc.

Is there any suggestions to improve the model like using different prior settings?

data {
  int<lower=1> N;
  int<lower=0,upper=1> y1[N];
  int<lower=1> K1; 
  matrix[N,K1] x1;
   vector<lower=0>[N] weights; 
}

parameters {
  real alpha1;
  vector[K1] beta;
   
}


model {
  
  vector[N] mu = alpha1+ x1 * beta;
  
  //priors
  beta ~ normal(0, 10);
  alpha1 ~ normal(0, 10);
   
  
     for (i in 1:N){
     target += weights[i] * bernoulli_logit_lpmf(y1[i] | mu[i]);;
     
       
     }
  
}

The warnings I got are as follows:

Thank you!!

Here are a couple ideas. First check that the weights are in the data{} block. I’d also suggest trying all weights equal to 1 to see if it works and compare to a “standard” unweighted logistic regression. If that all works, then make sure the weights you use sum to the sample size and not the population size. In other words, have the average weight value be 1. If you forgive the self-promotion, check out this package we are working on (if you are using R) GitHub - RyanHornby/csSampling

1 Like

Thank you I will look in to the package you suggested. I included the weights in my original code but forgot to include here.

Hi Ben,

I wanted to return to this comment and one you made similarly in a response to my post many moons ago. I tend to agree that using survey weights seems a violation of the likelihood principle. Similar arguments about violation of the likelihood principle have been leveled agains the use of the p-value, but in any case, it seems that Bayesians should be consistent about this. Having said that, are you aware of any references that explicitly raised the issue of the likelihood principle in the context of sampling weights - namely, as you say, one is conditioning on data that was never observed.

Thanks

David

I realize the question was directed at Ben, but I think the most elegant fully Bayesian approach is jointly modelling the selection probabilities and the outcome of interest. These recent papers show you can do it for glm type models, but you have to use custom stan code. I’m pretty sure the authors would be happy to share.

Luis G. León-Novelo. Terrance D. Savitsky. “Fully Bayesian estimation under informative sampling.” Electron. J. Statist. 13 (1) 1608 - 1645, 2019. Fully Bayesian estimation under informative sampling

2 Likes

@mrwilli: Do you know if the surhors address the issue of model feedback? AFAIK one cannot set up a fully bayesian model that jointly estimates weights and exposure outcome associations without bias.
See e.g. the discussion here: On Bayesian estimation of marginal structural models - PMC or here Model feedback in Bayesian propensity score estimation - PubMed

(I am ignoring the issue of weights not being Bayesians. There are some types of selection bias that MRP/standardization can’t deal with and that can only be dealt with with weights)

1 Like

Good question. I think it’s different for survey weights because they are observed, so this framework treats them as data and smooths them. Whereas for propensity weights they have to be estimated from binary indicators. So I have not seen this issue of model feedback for comodelling survey weights.

That being said we just worked on using a survey sample to come up with propensity weights for a convenience sample and I now wonder if we should check for this issue.

I agree that model feedback is not the problem if one tries to estimate the expected response in a population, like in opinion polls.
I was asking because the paper used the association between an exposure and an outcome as example. For such analyses model feedback os a problem.

Model feedback introduces bias because in a joint weight and outcome model, weight parameters that maximize the (weighted) likelihood of the outcome are preferred. As a result, parameters of the weight and outcome model will be biased.