Survey weights in brms/stan - Simulation based on design effect, feedback sought!

I’m aware of several discussions/threads on the forum discussing complex survey designs and what the weights are doing in brms (e.g., some insightful talk here: Are complex surveys feasible in brms?).

I have some projects in which I would like to get population-level estimates for certain parameters, and in frequentist approaches this can be achieved using survey weights and packages such as survey with the svyglm function. I am fully aware that survey weights are not considered ‘fully Bayesian’ and I will also be running assessments using full Bayesian MRP. However, there are some limits to what I can do with MRP, including correcting for biases in the sample for variables that I cannot combine into a workable post-stratification table. I also like to compute some follow-up parameters related to effect sizes that I cannot do with svyglm, nor with MRP, but could do with some form of ‘weighted regression’ using brms very nicely.

The weighting function in brm is essentially frequency weight, and therefore parameters estimated using survey weights (e.g., inverse probability weights) have accurate point estimates, but their error is overly optimistic as it isn’t corrected for the weighting procedure that would reduce the effective sample size.

Based on a comment from @simonbrauer I have experimented with ‘correcting’ the weights I feed into brms by calculating a general design effect, and then essentially penalising the sample size by dividing the weights by that design effect. I present some simulations below that show how this approach compares with svyglm (which I take to be ‘correct’ from a frequentist standpoint), and with typical weighting in brms (which does not behave as desired) - it seems that it works quite nicely in the sense that the errors using modified weights appropriately track what happens with proper survey weight analyses from svyglm, and I also show how the approach ‘correctly retrieves’ the true population value just as much as svyglm, and far more than typical brm weighting as the sample diverges from the population by increasing amounts.

Although I understand this is not ‘fully Bayesian’, I am wondering the extent to which it can be considered at least a reasonable approach. Note that I will often supplement these analyses with MRP estimates, so I won’t usually be just drawing straight conclusions only from these estimates in isolation - sort of like a multipronged approach to see how estimates converge using different methods.

Here is a plot showing how often the 95% intercept intervals contain the true parameter value across many simulations with different levels of bias (and hence weighting correction) in the data:

Here is a further plot showing how the typical widths of the 95% intervals with typical weighting are far too narrow, but the modified weighting approach closely tracks the width of the 95% confidence intervals from svyglm:

I’ve also uploaded all the code for these simulations, along with RData files so that you don’t actually have to spend the time running the sims, if you would like to check out how I have done this. GitHub - Jimbilben/Survey-Weighting-Simulation: This project looks at ways to calibrate brms weights to enable something similar to survey weights in a Bayesian or 'pseudo-Bayesian' regression model

I would really appreciate people’s thoughts on this. In particular I know there are lots of people who have been interested in or work in surveys at the moment e.g., @Corey_Sparks @Guido_Biele @bgoodri @jonah @lauren @maxbiostat @maurosc3ner (apologies if this is egregious tagging but your comments were all informative in previous posts about weighting/MRP!).

The key goal with this is really to have access to the full Bayesian suite of posterior predictions, adding priors and so on, but not generating estimates that are excessively optimisitic with regards to the precision we can reach, by incorporating some penalties based upon more typical survey weighting procedures - from my simulations this does seem to work.

4 Likes

Here are two further plots showing the approach also seems to work for a bernoulli/binomial model. In this case, with highly unrepresentative data it seems like it starts to perform slightly better than svyglm. This might be
because by adding a prior I was able to constrain what values would be entertained by the model - keeping them more realistic and also meaning that the width of the posterior only widens to a certain degree while still managing to contain the true value.


1 Like

We put out a simple package that can fit a weighted model in brms and postprocess using the gradient functions via rstan.

It’s just r code wrappers for rstan and the survey package. There’s also a wrapper for brms. It needs some more work and I hope to update it soon with some more options for estimating the sandwich matrices. Please reach out if you have any issues. For a simple glm it should work fine but we’ve only had a handful of testers. There are a few examples in the bundled rhelp files. The key is that we need rstan to access the gradients via autodiff.

5 Likes

Thanks @mrwilli that’s great news. The brms wrapper threw me an error that it would not continue because the mean of the weights does not equal 1. The mean is like 1.000000000045 or something like that. Is there a way to make the model just run with this ‘essentially is 1’ mean of the weights?

If you would be able to see what you think is going wrong here it would be realy appreciated. In the first instance, the model does not run because there is a miniscule remainder in the mean of the weights not exactly equalling zero.

But if I just set all the weights to 1 to make sure the mean is 1, I get a different error in any case:

Error in object@.MISC$stan_fit_instance$unconstrain_pars(pars) :
mismatch in number dimensions declared and found in context; processing stage=parameter initialization; variable name=b; dims declared=(1); dims found=()

The code I’m using is as follows, and I’m attaching the test_data example:

library(csSampling)
library(rstan)
library(brms)
library(survey)
rstan_options(auto_write = TRUE)

svy_weight_des <- svydesign(ids = ~1, data = test_data, weights = test_data$weight)

# just looking to make a simple model estimating the weighted average
svy_brms <- cs_sampling_brms(svydes = svy_weight_des, 
                             brmsmod = brmsformula(score | weights(weight) ~ 1, center = FALSE), 
                             data = test_data,
                             family = gaussian())

mean(test_data$weight)
mean(test_data$weight) == 1
mean(test_data$weight) %% 1 # there is a tiny remainder that means it != 1

# make sure the weights average to 1!
test_data$weights_check <- 1

svy_weight_des2 <- svydesign(ids = ~1, data = test_data, weights = test_data$weights_check)

svy_brms2 <- cs_sampling_brms(svydes = svy_weight_des2, 
                             brmsmod = brmsformula(score | weights(weights_check) ~ 1, center = FALSE), 
                             data = test_data,
                             family = gaussian())

# Throws a different error:
# Error in object@.MISC$stan_fit_instance$unconstrain_pars(pars) : 
# mismatch in number dimensions declared and found in context;
# processing stage=parameter initialization; variable name=b; dims declared=(1); dims found=()

example data.RData (23.4 KB)

Any help would be really appreciated!

Thanks for your patience and sharing your example!

The issue is a helper function that coerces the parameter draws into something that can be input in the rstan::unconstrain_pars function. (This has been the hardest part of the code for me!) The current helper drops dimensions of 1, but needs to be fixed. I think if you have at least one predictor in addition to an intercept - beta vector is dim 2, it should work.

For the mean != 1 problem try this

test_data$weights_check <- test_data$weight - (sum(test_data$weight) - length(test_data$weight))/length(test_data$weight)

The goal is to get the sum of the weights equal to the length or sample size. usually I do a ratio adjustment, but this difference is so small, it seems to only work as an additive adjustment.

Using this model with more predictors, I got it to run for me:

svy_brms2 <- cs_sampling_brms(svydes = svy_weight_des2, 
                             brmsmod = brmsformula(score | weights(weights_check) ~ age + sex, center = FALSE), 
                             data = test_data,
                             family = gaussian())

plot(svy_brms2)

Below is the current helper function and then the new one that needs to be swapped out. I’ll open an issue on github to motivate this fix.

list_2D_row_subset

function (nmlist, rindex) 
{
    temp_list <- list()
    for (k in 1:length(nmlist)) {
        ldim <- length(dim(nmlist[[k]]))
        lcommas <- paste(rep(",", ldim - 1), collapse = " ")
        eval(parse(text = paste("temp_list$", names(nmlist)[k], 
            " <- ", "(nmlist$", names(nmlist)[k], 
            ")[rindex", lcommas, "]", sep = "")))
    }
    return(temp_list)
}

We need to replace it with the following which is more careful for arrays with size dimensions of 1

function (nmlist, rindex) 
{
    temp_list <- list()
    for (k in 1:length(nmlist)) {
	  tmpdim <- dim(nmlist[[k]])
        ldim <- length(tmpdim)
        lcommas <- paste(rep(",", ldim - 1), collapse = " ")
	 #copy over to new list - drop = FALSE retains ALL dimensions
        eval(parse(text = paste("temp_list$", names(nmlist)[k], 
            " <- ", "(nmlist$", names(nmlist)[k], 
            ")[rindex", lcommas, ",drop = FALSE]", sep = "")))
       #drop only first dimension of array - not others of size 1
	 if(ldim > 1){
	  eval(parse(text = paste("temp_list$", names(nmlist)[k], 
            " <- ", "array(temp_list$", names(nmlist)[k], ", dim = tmpdim[-1])", sep = "")))
       }	  	
    }
    return(temp_list)
}

1 Like

You can ignore the message above - I looked in the github thing and you wrote there it is intended to circumvent the intercept only problem. Thank you @mrwilli I’ll hopefully be able to test it out later today

Do you have any thoughts on how I could modify the code I pasted above to incorporate this new helper function and get it to run? I tried just defining the function in my R environment but it didn’t help - which is probably to be expected! But not sure where/how I should define it to make it work with the other functions.

I’m showing my newness to making packages, but yes we can’t simply rename an associated new helper function in our local session and have the cs_sampling function recognize it. I forget the term, but csSampling namespace is locked or protected, so we’d have to update and reinstall the package. As an alternative, try the code below with different names. It’s bit of a patch for now, just rewriting the functions, but it should work.

list_2D_row_subset_DROP <- function (nmlist, rindex) 
{
    temp_list <- list()
    for (k in 1:length(nmlist)) {
	  tmpdim <- dim(nmlist[[k]])
        ldim <- length(tmpdim)
        lcommas <- paste(rep(",", ldim - 1), collapse = " ")
	 #copy over to new list - drop = FALSE retains ALL dimensions
        eval(parse(text = paste("temp_list$", names(nmlist)[k], 
            " <- ", "(nmlist$", names(nmlist)[k], 
            ")[rindex", lcommas, ",drop = FALSE]", sep = "")))
       #drop only first dimension of array - not others of size 1
	 if(ldim > 1){
	  eval(parse(text = paste("temp_list$", names(nmlist)[k], 
            " <- ", "array(temp_list$", names(nmlist)[k], ", dim = tmpdim[-1])", sep = "")))
       }
      #if only had 1 dim which is the MCMC draw, make a double (no dim), rather than an array of dim 1 or 0
	if(ldim == 1){
	  eval(parse(text = paste("temp_list$", names(nmlist)[k], 
            " <- ", "as.double(temp_list$", names(nmlist)[k], ")", sep = "")))
	}	  	
    }
    return(temp_list)
}


cs_sampling_DROP <- function (svydes, mod_stan, par_stan = NA, data_stan, ctrl_stan = list(chains = 1, 
    iter = 2000, warmup = 1000, thin = 1), rep_design = FALSE, 
    ctrl_rep = list(replicates = 100, type = "mrbbootstrap"), 
    sampling_args = list()) 
{
    require(rstan)
    require(survey)
    require(plyr)
    require(pkgcond)
    if (rep_design) {
        svyweights <- svydes$pweights
    }else {
        svyweights <- weights(svydes)
    }
    if (is.null(svyweights)) {
        if (!is.null(weights(data_stan))) {
            stop("No survey weights")
        }
    }
    if (is.null(weights(data_stan))) {
        if (!is.null(svyweights)) {
            warning("No stan data weights, using survey weights instead")
            data_stan$weights = weights(svydes)
        }
    }
    if (!isTRUE(all.equal(as.numeric(weights(data_stan)), as.numeric(svyweights)))) {
        stop("Survey weights and stan data weights do not match")
    }
    if (mean(weights(data_stan)) != 1) {
        stop("Mean of the weights is not 1")
    }
    print("stan fitting")
    out_stan <- do.call(sampling, c(list(object = mod_stan, data = data_stan, 
        pars = par_stan, chains = ctrl_stan$chains, iter = ctrl_stan$iter, 
        warmup = ctrl_stan$warmup, thin = ctrl_stan$thin), sampling_args))
    par_samps_list <- rstan::extract(out_stan, permuted = TRUE)
    if (anyNA(par_stan)) {
        par_stan <- names(par_samps_list)[-length(names(par_samps_list))]
    }
    par_samps <- as.matrix(out_stan, pars = par_stan)
    for (i in 1:dim(par_samps)[1]) {
        if (i == 1) {
            upar_samps <- unconstrain_pars(out_stan, list_2D_row_subset_DROP(par_samps_list, 
                i))
        }
        else {
            upar_samps <- rbind(upar_samps, unconstrain_pars(out_stan, 
                list_2D_row_subset_DROP(par_samps_list, i)))
        }
    }
    row.names(upar_samps) <- 1:dim(par_samps)[1]
    upar_hat <- colMeans(upar_samps)
    Hhat <- -1 * optimHess(upar_hat, gr = function(x) {
        grad_log_prob(out_stan, x)
    })
    if (rep_design == TRUE) {
        svyrep <- svydes
    }
    else {
        svyrep <- as.svrepdesign(design = svydes, type = ctrl_rep$type, 
            replicates = ctrl_rep$replicates)
    }
    print("gradient evaluation")
    rep_tmp <- withReplicates(design = svyrep, theta = grad_par, 
        stanmod = mod_stan, standata = data_stan, par_hat = upar_hat)
    Jhat <- vcov(rep_tmp)
    Hi <- solve(Hhat)
    V1 <- Hi %*% Jhat %*% Hi
    R1 <- chol(V1, pivot = TRUE)
    pivot <- attr(R1, "pivot")
    R1 <- R1[, order(pivot)]
    R2 <- chol(Hi, pivot = TRUE)
    pivot2 <- attr(R2, "pivot")
    R2 <- R2[, order(pivot2)]
    R2i <- solve(R2)
    R2iR1 <- R2i %*% R1
    upar_adj <- aaply(upar_samps, 1, DEadj, par_hat = upar_hat, 
        R2R1 = R2iR1, .drop = TRUE)
    for (i in 1:dim(upar_adj)[1]) {
        if (i == 1) {
            par_adj <- unlist(constrain_pars(out_stan, upar_adj[i, 
                ])[par_stan])
        }
        else {
            par_adj <- rbind(par_adj, unlist(constrain_pars(out_stan, 
                upar_adj[i, ])[par_stan]))
        }
    }
    row.names(par_adj) <- 1:dim(par_samps)[1]
    colnames(par_samps) <- colnames(par_adj)
    rtn = list(stan_fit = out_stan, sampled_parms = par_samps, 
        adjusted_parms = par_adj)
    class(rtn) = c("cs_sampling", class(rtn))
    return(rtn)
}

cs_sampling_brms_DROP <- function (svydes, brmsmod, data, family, par_brms = NA, prior = NULL, 
    stanvars = NULL, knots = NULL, ctrl_stan = list(chains = 1, 
        iter = 2000, warmup = 1000, thin = 1), rep_design = FALSE, 
    ctrl_rep = list(replicates = 100, type = "mrbbootstrap"), 
    stancode_args = list(), standata_args = list(), sampling_args = list()) 
{
    stancode <- do.call(make_stancode, c(list(brmsmod, data = data, 
        family = family, prior = prior, stanvars = stanvars, 
        knots = knots), stancode_args))
    print("compiling stan model")
    mod_brms <- stan_model(model_code = stancode)
    data_brms <- do.call(make_standata, c(list(brmsmod, data = data, 
        family = family, prior = prior, stanvars = stanvars, 
        knots = knots), standata_args))
    return(cs_sampling_DROP(svydes = svydes, mod_stan = mod_brms, 
        par_stan = par_brms, data_stan = data_brms, rep_design = rep_design, 
        ctrl_rep = ctrl_rep, ctrl_stan = ctrl_stan, sampling_args))
}

These two examples now worked for me

svy_brms2_DROP <- cs_sampling_brms_DROP(svydes = svy_weight_des2, 
                             brmsmod = brmsformula(score | weights(weights_check) ~ age + sex, center = FALSE), 
                             data = test_data,
                             family = gaussian())

plot(svy_brms2_DROP)
svy_brms3_DROP <- cs_sampling_brms_DROP(svydes = svy_weight_des2, 
                             brmsmod = brmsformula(score | weights(weights_check) ~ 1, center = FALSE), 
                             data = test_data,
                             family = gaussian())

plot(svy_brms3_DROP)

In the previous post, the code I suggested worked for the regression coefficients but not sigma. Hence the extra condition in list_2D_row_subset_DROP for ldim == 1

2 Likes

Nice one! Thank you yes this does run!

No problem that this isn’t fully updated in the package yet - I can probably just bundle it as a package on my local machine so it is easier to call, until it is fully in your package.

Do you know how far the support for different types of models extends with this? For example, would it appropriately calculate parameters for a repeated measures design where you weight each subject according to population representativeness, but each subjects appears twice e.g., before and after some intervention?

Good question. One approach is to fit an independence model repeating the same individuals over time and making sure each is represented in the same cluster ID in the survey design. This approach would estimate a “robust” sandwich covariance matrix. It treats the repeated measures as a nuisance correlation, but the regression coefficients and other parameters have the variance adjusted. If on the other hand, we’d want to do a variance decomposition and treat the individual as a group level random effect, this gets trickier and we’d probably want to weight the random effect distribution for the individual, leading to multiple weights. I don’t think it can be directly done in brms and adjusted. We tried the estimation part in a couple papers and have recently tried the adjustment for a simple one-way ANOVA. The cs_sampling adjustment didn’t work “out of the box” and had to be tweaked. If you have a motivating example, we can try to work through it!

I would mostly actually be thinking of quite simple cases: we recruit a sample and have weights for each respondent. The repeated measures component might just be that we ask the same person multiple questions, and so we want to fit a random intercept for participant, as well as an effect of - for example - question framing.

I’ve confirmed that with my relatively simple approach above of just dividing by the design effect, it can correctly get at repeated measures effects but it would be nice to be able to do this with another package.

@mrwilli did you have any thoughts about using svybrms in a repeated measures fashion? I am imgainging a case where we recruit people all at once and ask them how much they support some different policies. In an unweighted version this would use something like the simple formula:

support ~ policy + (1 | respondent)

So it is repeated measures, but everyone is collected at once and has a survey weight for recruitment at that time - it doesn’t need to be corrected e.g., for being sampled again at a later date or something like that. Can svybrms (or just survey) handle these quite simple repeated measures designs?