Site icon R-bloggers

What Do Tunisians Really Think About President Kais Saied?

[This article was first published on R on Robert Kubinec, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This Rmarkdown document contains code and data for the Kais Saied survey experiment by Robert Kubinec and Amr Yakout. This survey was fielded using an online panel in Tunisia during August of 2023. The survey experiment involved providing half the sample a direct question about whether they supported Kais Saied’s move to suspend parliament and centralize power in his hands, and the other half of the sample had a randomized response question designed to allow them to answer truthfully without being identified. This experiment was pre-registered and the pre-registration can be accessed from this link.

This document was drafted to allow people to verify our results using the raw data. The survey data included has been stripped of any identifying information about respondents, but it does include the actual answers from the survey so that all analyses can be reproduced. You can access the code file and the raw data from this Github repository: https://github.com/saudiwin/saudiwin.github.io/tree/sources/content/post (file is entitled kais_saied_results.Rmd and the data file is in the data/ subfolder). The survey data can be found in the data/ folder as kais_saied_survey.csv.

This Quarto document includes embedded R code that loads the survey data and estimates the Stan model for proportion of people who support Kais Saied, and compares it to the direct question to see if there is sensitive survey bias. At present there are a total of 835 completed responses.

All questions about the analysis should be directed to Robert Kubinec at rmk7@nyu.edu.

Robustness check: time of completion

In this section I show some general statistics for the data. In particular, I look at how long it takes people to complete the survey. In the plot below I show the average duration in minutes over time for each survey response. The plot shows that most people took the survey in about 10 to 30 minutes:

survey_data %>% 
  filter(`Duration (in seconds)`<10000) %>% 
  mutate(Duration=`Duration (in seconds)`/60) %>% 
  ggplot(aes(y=Duration,x=StartDate)) +
  geom_point() +
  scale_y_log10() +
  geom_smooth()

The plot also shows that data collection started in late July and continued until mid-August. Data came in fairly regularly on a daily basis.

Survey Design

In the survey we implemented a form of a randomized response question that is designed to help people report answers on a survey that they may be embarrassed about or in danger if they report. The question was originally developed for private medical information like abortion or HIV infections, and has been extended to cover crime, drug use and corruption reporting. See Hout and Heijden (2002), Gingerich et al. (2016) and Blair, Imai, and Zhou (2015) for more information.

The theory behind a randomized response is to inject randomness into the survey question. Theoretically it is similar to approaches in cryptography that use random numbers to encode online data using keys. In our case, we used naturally ocurring randomness related to people’s birthdays. Because we did not ask for respondents’ birth dates, and the month in which one is born is essentially random, we can use that naturally occurring randomness to encode their responses.

To test for how much we could encourage people to respond truthfully, we randomly assigned each respondent to receive one of the two questions below:

Direct Question

  1. Do you oppose President Kais Saied’s moves to change Tunisia’s constitution and close the parliament?
    1. Yes
    2. No

Randomized Response

We understand that politics in Tunisia is sensitive right now. This question is worded so that you can tell us what you think but still protect your privacy. Because we don’t know when your mother was born, we also won’t know for sure your political opinion.

  1. My mother’s birth date is in January, February or March.
  2. I oppose President Kais Saied’s moves to change Tunisia’s constitution and close the parliament.

Please pick the answer that best represents whether these statements are true of you:

  1. Both statements are true OR neither is true.
  2. One of the two statements is true.

While explaining how this question works is beyond the scope of this note, we again refer to the linked paper above for the statistical mechanics. As long as the respondent reads and follows the instructions, they can report their opposition (or support) for Kais Saied’s power grab without us being able to know their true response. Essentially, we can’t separate their answer from whether their mother’s birthday is in a given month, and as a result, the individual answers are encoded. For any respondent, we have no idea whether they oppose Saied or not.

However, the beauty of this method is that even though we don’t know any one person’s response, we can estimate how many in people in total support or oppose Saied. When we aggregate responses, we can take into account that on average \(\frac{1}{4}\) of our respondents will have mothers who were born in those months. Using statistical models that we estimate below, we can find out how many people truly oppose Saied, and also how many people appear to be under-reporting their opposition relative to the direct question.

Sensitivity Estimates

In this section we show several different ways that we calculate the true number of people who are opposed to Kais Saied. We first use a simple Bayesian calculation based on the beta distribution that was specified in the pre-registration linked above. We then look at other existing R packages that are used for randomized response questions and then we implement a custom Bayesian model that allows us to jointly model both the randomized response and the direct question. We also use this custom model to allow us to adjust for biases in using online panels that may not be fully representative.

# proportion selecting birthday response

lambda <- table(survey_data$kais_rr)
lambda <- lambda[2]/sum(lambda)

N <- sum(as.numeric(!is.na(survey_data$kais_rr) & survey_data$FL_63_DO_saied_sensitive==1))

to_stan_data <- list(success_cm=N*((lambda - .75)/(-.5)),
                     fail_cm= N*(1-((lambda - .75)/-.5)),
                     success_d=sum(as.numeric(survey_data$kais_direct=="Yes" & survey_data$FL_63_DO_saied_nonsensitive==1),na.rm=T),
                     fail_d=sum(as.numeric(survey_data$kais_direct=="No" & survey_data$FL_63_DO_saied_nonsensitive==1),na.rm=T))

fit_mod <- sen_mod$sample(chains=4, cores=4, data=to_stan_data,refresh=0)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
fit_mod$summary()
## # A tibble: 5 × 10
##   variable           mean   median     sd    mad       q5      q95  rhat ess_b…¹
##   <chr>             <num>    <num>  <num>  <num>    <num>    <num> <num>   <num>
## 1 lp__           -790.    -790.    1.21   1.00   -792.    -789.     1.00   2093.
## 2 pi_hat            0.512    0.512 0.0248 0.0244    0.471    0.553  1.00   4310.
## 3 d_hat             0.316    0.316 0.0227 0.0227    0.279    0.353  1.00   4612.
## 4 d_hat_binomial    0.315    0.315 0.0233 0.0238    0.277    0.353  1.00   4161.
## 5 estimand          0.196    0.196 0.0332 0.0337    0.142    0.251  1.00   4350.
## # … with 1 more variable: ess_tail <num>, and abbreviated variable name
## #   ¹​ess_bulk

We can plot the estimated bias as the difference between the average proportion of those responding affirmative to the direct question and the calculated level of opposition to Saied from the randomized response question:

estimand <- c(fit_mod$draws("estimand"))

mcmc_dens(fit_mod$draws("estimand")) + 
  theme_tufte()

Our basic model shows that approximately 19.6% of respondents reported anti-Saied preferences in the randomized response method versus the direct question. The 5% – 95% uncertainty interval is from 14.2 to 25.1 .

The full posterior density of the bias-corrected estimate is:

mcmc_dens(fit_mod$draws("pi_hat")) + theme_tufte()

Comparing the two estimates as intervals (pi_hat is the true level of opposition to Saied and d_hat is the proportion from the direct question):

mcmc_intervals(fit_mod$draws(c("pi_hat","d_hat"))) + theme_tufte()

Alternative Estimation: RRreg

In this section we use the R package RRreg to check these calculations using a more traditional means for estimating sensitive proportions.

est_rreg <- RRuni(response=as.numeric(comp_data_rr$kais_rr=="Both statements are true OR neither is true." & comp_data_rr$FL_63_DO_saied_sensitive==1),
                  p=.25,
                  model="Crosswise",MLest = TRUE)

summary(est_rreg)
## Crosswise Model with p = 0.25
## Sample size: 391
## 
##    Estimate   StdErr      z  Pr(>|z|)    
## pi 0.512788 0.050633 10.128 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
s1 <- summary(est_rreg)

Let’s compare these estimates with plots:

est_beta <- fit_mod$summary()

est_data <- tibble(estimate=c(s1$coefficients[,1],
                              est_beta$median[est_beta$variable=="pi_hat"]),
                   low=c(s1$coefficients[,1]-1.96*(s1$coefficients[,2]),
                         est_beta$q5[est_beta$variable=="pi_hat"]),
                   high=c(s1$coefficients[,1]+1.96*(s1$coefficients[,2]),
                         est_beta$q95[est_beta$variable=="pi_hat"]),
                   estimator=c("RRreg","Bayes"))

est_data %>% 
  ggplot(aes(y=estimate,x=estimator)) +
  geom_linerange(aes(ymin=low,
                     ymax=high)) +
  geom_point() +
   theme_tufte()

We can see that they are similar estimates, although the R package has higher uncertainty than the simple Bayesian estimate.

Alternative Estimator: rr

We can also use the R package rr which uses the EM algorithm as opposed to RRreg which uses more conventional optimization.

comp_data_rr$gender_two <- factor(comp_data_rr$gender,
                                  exclude="Other:")
comp_data_rr <- filter(comp_data_rr, !is.na(gender_two))
# need to switch the outcome
comp_data_rr$kais_rr_num_switch <- 1 - comp_data_rr$kais_rr_num
rr_est <- rrreg(kais_rr_num_switch ~ as.character(gender_two), 
                data=comp_data_rr,
                      p=.75,
                      design='mirrored')

# predict sensitive trait

pred_rr <- predict(rr_est,quasi.bayes = T)
mean(pred_rr$est)
## [1] 0.5083205
mean(pred_rr$est[comp_data_rr$gender_two=="Male"],na.rm=T)
## [1] 0.5125389
mean(pred_rr$est[comp_data_rr$gender_two=="Female"],na.rm=T)
## [1] 0.5018136
pred_rr_mean <- predict(rr_est,quasi.bayes = T,avg=T)

This number seems very close to our other estimators. Again, the uncertainty seems larger than the simple Bayesian method. Minimal gender differences over who is more or less likely to oppose Saied.

We can plot this estimate against the others:

est_data <- tibble(estimate=c(s1$coefficients[,1],
                              est_beta$median[est_beta$variable=="pi_hat"],
                              pred_rr_mean$est),
                   low=c(s1$coefficients[,1]-1.96*(s1$coefficients[,2]),
                         est_beta$q5[est_beta$variable=="pi_hat"],
                         pred_rr_mean$ci.lower),
                   high=c(s1$coefficients[,1]+1.96*(s1$coefficients[,2]),
                         est_beta$q95[est_beta$variable=="pi_hat"],
                         pred_rr_mean$ci.upper),
                   estimator=c("RRreg","Bayes","EM"))

est_data %>% 
  ggplot(aes(y=estimate,x=estimator)) +
  geom_linerange(aes(ymin=low,
                     ymax=high)) +
  geom_point() + theme_tufte()

BRMS Model

Finally, we define a brms custom family for the randomized response model using the confusion matrix approach of Hout and Heijden (2002). This model is a variety of the Bernoulli distribution that takes into account the known probabilities of obtaining the true response. It jointly models both the randomized response model and the direct question, allowing us to estimate the bias directly rather than post-estimation. By implementing it in brms, we can make use of brms features like multilevel regression to allow us to do post-stratification adjustment of the survey responses with population data from Tunisia’s 2014 census. This is our preferred specification.

library(brms)

stan_funs <- '
  real sens_reg_lpmf(int y, real mu, real bias, matrix P, int T) {
  
    // generalized RR model from van der Hout and van der Heijden (2002)
    // also see R package RRreg
    
    real out;
    // need to impose a constraint on bias where it cannot be larger than mu
    real bias_trans = mu * bias;
    
    if(T==1) {
    
      // treatment distribution (crosswise model)
      
      if(y==1) {
      
        out = P[2,1] * (1 - mu) + P[2,2] * mu;
      
      } else if(y==0) {
      
       out = P[1,1]*(1-mu) + P[1,2]*mu;
      
      }

    } else if (T==0) {
    
      // control = direct question
    
      if(y==1) {
      
        out = mu - bias_trans;
      
      } else if(y==0) {
      
        out = (1 - mu) + bias_trans;
      
      }
    
    }
    
    return log(out); 

  }
  
  int sens_reg_rng(real mu, real bias, matrix P, int T) {
  
    real bias_trans = mu*bias;
  
    if(T==1) {
    
      return bernoulli_rng(P[2,1] * (1 - mu) + P[2,2] * mu);
      
    } else {
    
      return bernoulli_rng(mu - bias_trans);
    
    }

  }'

# define custom family

family_sens_reg <- custom_family("sens_reg",
                                 dpars=c("mu","bias"),
                                 links=c("logit","logit"),
                                 type="int",
                                 lb=c(NA,NA),
                                 ub=c(NA,NA),
                                 vars=c("P","vint1[n]"))

# define log-likelihood

log_lik_sens_reg <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  y <- prep$data$Y[i]
  treatment <- prep$data$vint1[i]
  bias <- brms::get_dpar(prep, "bias", i = i)
  
  bias_trans <- bias*mu
  
  if(treatment==1) {
    
    if(y==1) {
    
      return(log(P[2,1] * (1 - mu) + P[2,2] * mu))
  
    } else {
      
      return(log(P[1,1]*(1-mu) + P[1,2]*mu))
    
    }
    
  } else {
    
    if(y==1) {
      
      return(log(mu - bias_trans))
      
    } else {
      
      return(log((1 - mu) + bias_trans))
      
    }
    
  }
  

}

# define posterior predictions

posterior_predict_sens_reg <- function(i, prep, ...) {
  
  mu <- brms::get_dpar(prep, "mu", i = i)
  bias <- brms::get_dpar(prep, "bias", i = i)
  y <- prep$data$Y[i]
  treatment <- prep$data$vint1[i]
  
  bias_trans <- mu*bias

  if(treatment==1) {
    
    out <- rbinom(n=length(mu),size=1,prob=P[2,1] * (1 - mu) + P[2,2] * mu)
    
  } else {
    
    out <- rbinom(n=length(mu),size=1,prob=mu - bias_trans)
    
  }
  
  return(out)
  
}

# define posterior expectation (equal to latent variable pi)

posterior_epred_sens_reg <- function(prep,...) {

  mu <- brms::get_dpar(prep, "mu")
  bias <- brms::get_dpar(prep, "bias")

  mu

}

posterior_epred_bias_sens_reg <- function(prep,...) {

  mu <- brms::get_dpar(prep, "mu")
  bias <- brms::get_dpar(prep, "bias")

  bias*mu

}

P <- getPW("Warner",p=.25)

all_stanvars <- stanvar(x=P,block = "data") + 
  stanvar(scode=stan_funs,block="functions")

survey_data$age_cat_order <- ordered(survey_data$age_cat)

fit1 <- brm(bf(kais_combined | vint(treatment) ~ gender*mo(age_cat_order) + (1|gov)), 
            data=survey_data,
            family=family_sens_reg,
            stanvars=all_stanvars,
            prior=prior(beta(1,1),class="bias") + 
              prior(normal(0,5), class="b"),
            chains=2,cores=2,control=list(max_treedepth=11,
                                          adapt_delta=0.95),
            backend = "cmdstanr")

pp_check(fit1, type="bars",ndraws=500)

loo(fit1)
## 
## Computed from 2000 by 804 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -534.4  7.8
## p_loo        10.9  0.4
## looic      1068.9 15.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     803   99.9%   702       
##  (0.5, 0.7]   (ok)         1    0.1%   1529      
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
gov_cats <- ranef(fit1,groups = "gov")$gov[,,1] %>% as_tibble %>% 
  mutate(level=row.names(ranef(fit1,groups = "gov")$gov[,,1]))

Plot this additional type of estimation:

get_latent_draws <- posterior_epred(fit1)

get_latent_est <- tibble(median=median(get_latent_draws),
                         high=quantile(apply(get_latent_draws, 1, median),.95),
                         low=quantile(apply(get_latent_draws, 1, median),.05))

# get estimate of bias

prep_bias <- prepare_predictions(fit1)

bias_trans <- posterior_epred_bias_sens_reg(prep_bias)

bias_trans_est <- tibble(bias_est_low=quantile(apply(bias_trans,1,median),.05),
                         bias_est=median(apply(bias_trans,1,median)),
                         bias_est_high=quantile(apply(bias_trans,1,median),.95))

est_data <- tibble(estimate=c(s1$coefficients[,1],
                              est_beta$median[est_beta$variable=="pi_hat"],
                              pred_rr_mean$est,
                              get_latent_est$median),
                   low=c(s1$coefficients[,1]-1.96*(s1$coefficients[,2]),
                         est_beta$q5[est_beta$variable=="pi_hat"],
                         pred_rr_mean$ci.lower,
                         get_latent_est$low),
                   high=c(s1$coefficients[,1]+1.96*(s1$coefficients[,2]),
                         est_beta$q95[est_beta$variable=="pi_hat"],
                         pred_rr_mean$ci.upper,
                         get_latent_est$high),
                   estimator=c("RRreg","Bayes","EM","Bayes_Reg"))

est_data %>% 
  ggplot(aes(y=estimate,x=estimator)) +
  geom_linerange(aes(ymin=low,
                     ymax=high)) +
  geom_point() + theme_tufte()

knitr::kable(bias_trans_est)
bias_est_low bias_est bias_est_high
0.1155721 0.2164725 0.3187978

Again, we see that all of the estimates are quite similar. In the last section we compare them directly with a simulation, which shows that the two Bayesian estimators are both accurate and have good coverage, while RRreg is accurate but understates uncertainty in general.

Adjustment with MRP

Given that the data come from an online panel, we want to use Tunisian 2014 census data to adjust these findings by age, sex and governorate. This will account for a considerable (though not all) amount of sample selection bias due to the online survey frame. We will use random effects to model the influence of the sample frame on the outcome.

head(census)
## # A tibble: 6 × 6
##   gender gov     pop age_cat age_cat_order    prop
##   <chr>  <chr> <dbl> <fct>   <ord>           <dbl>
## 1 Male   Tunis 17446 15-19   15-19         0.00120
## 2 Female Tunis 16984 15-19   15-19         0.00117
## 3 Male   Tunis 54142 20-24   20-24         0.00373
## 4 Female Tunis 51111 20-24   20-24         0.00352
## 5 Male   Tunis 48773 25-29   25-29         0.00336
## 6 Female Tunis 45668 25-29   25-29         0.00315

We will predict the sensitive trait for each census category, then stratify by summing over the proportion of population in each cell from the census data. We then plot this adjusted estimate.

# do for each condition, then average

census_treat <- mutate(census, treatment=1)
census_notreat <- mutate(census, treatment=0)

pred_data <- bind_rows(census_treat,
                       census_notreat) %>% 
  mutate(prop=prop/2)

tunisia_pred <- posterior_epred(fit1, newdata=pred_data) %>% 
  as_tibble %>% 
  mutate(draw=1:n()) %>% 
  gather(key="key",
         value="estimate",-draw) %>% 
  mutate(draw=paste0("V",draw)) %>% 
  spread(key="draw",value="estimate") %>% 
  mutate(key=as.numeric(stringr::str_remove(key, "V")))

tunisia_pred <- left_join(tunisia_pred, 
                          mutate(pred_data, key=1:n())) %>% 
  left_join(census) %>% 
  gather(key = "draw",
         value= "estimate",
         matches("V",ignore.case=F))

# aggregate to highest level

agg_est <- tunisia_pred %>% 
  group_by(draw) %>% 
  summarize(est_adj = sum(estimate * prop))

agg_est %>% 
  ggplot(aes(x=est_adj)) +
  geom_histogram() + theme_tufte()

quantile(agg_est$est_adj, c(0.05,.5,.95))
##        5%       50%       95% 
## 0.4198164 0.5116073 0.6027915

A difference of approximately 1 percent between the adjusted estimate and the naive number:

quantile(apply(get_latent_draws, 1, mean),c(0.05,.5,.95))
##        5%       50%       95% 
## 0.4337994 0.5237461 0.6165752

Next we can aggregate these results to look at gender, age and governorate in terms of relative levels of Saied support. Given that there is substantial uncertainty due to the sensitive question design, these results are somewhat noisy and should be interpreted with caution.

Gender

When we aggregate to the level of gender, we see that men are approximately 5% more likely than women to report opposition President Kais Saied.

# merge in census data
tunisia_pred %>% 
  group_by(gender,draw) %>% 
  mutate(prop_gender=pop/sum(pop)) %>% 
  summarize(est_gender=sum(estimate*prop_gender)) %>% 
  group_by(gender) %>% 
  summarize(median_gender=median(est_gender),
            low_gender=quantile(est_gender,.05),
            high_gender=quantile(est_gender,.95)) %>% 
  ggplot(aes(y=median_gender,
             x=gender)) +
  geom_pointrange(aes(ymin=low_gender,
                      ymax=high_gender)) +
  scale_y_continuous(labels=scales::percent) +
    labs("% Opposing Saied") +
  theme_tufte()

Age

The plot below shows aggregated opposition by age. Approximately 60 percent of the youngest age category (18 to 19 year olds) oppose President Saied while only 40 percent of the oldest age category (80 and over) oppose President Saied. This pattern is noticeably stronger than that for gender, though there is still considerable uncertainty.

tunisia_pred %>% 
  group_by(age_cat,draw) %>% 
  mutate(prop_age_cat=pop/sum(pop)) %>% 
  summarize(est_age_cat=sum(estimate*prop_age_cat)) %>% 
  group_by(age_cat) %>% 
  summarize(median_age_cat=median(est_age_cat),
            low_age_cat=quantile(est_age_cat,.05),
            high_age_cat=quantile(est_age_cat,.95)) %>% 
  ggplot(aes(y=median_age_cat,
             x=age_cat)) +
  geom_pointrange(aes(ymin=low_age_cat,
                      ymax=high_age_cat)) +
    scale_y_continuous(labels=scales::percent) +
  labs("% Opposing Saied") +
   theme_tufte() +
  theme(axis.text.x = element_text(angle=90))

Region

Finally, the plot below shows predicted opposition by region. In general, more rural districts like Kebili, Beja and Kairouan report more opposition, but these differences are quite noisy. Differences between regions are not especially pronounced in the data.

tunisia_pred %>% 
  group_by(gov,draw) %>% 
  mutate(prop_gov=pop/sum(pop)) %>% 
  summarize(est_gov=sum(estimate*prop_gov)) %>% 
  group_by(gov) %>% 
  summarize(median_gov=median(est_gov),
            low_gov=quantile(est_gov,.05),
            high_gov=quantile(est_gov,.95)) %>% 
  ggplot(aes(y=median_gov,
             x=reorder(gov,median_gov))) +
  coord_flip() +
    labs("% Opposing Saied") +
  geom_pointrange(aes(ymin=low_gov,
                      ymax=high_gov)) + theme_tufte()

Simulation Comparison: Univariate Vs. Regression

In this simulation we test the coverage and unbiasedness of the estimators, including our preferred Bayesian specification. The simulation shows that the estimators perform quite well at recovering the sensitive trait, and our Bayesian models have excellent coverage (the uncertainty they report is reasonable). We’ll use the regression spec as the DGP along with our observed parameters for this study.

library(parallel)

# confusion matrix

P <- getPW("Warner",.25)

# assume N = 500, small sensitivity of 0.05

theta <- 0.3
N <- 800
bias <- 0.1

sims <- 500

if(run_sim) {
  
  over_sims <- lapply(1:sims, function(s) {
  
  # assign half to treatment, half to control
  
  treatment <- as.numeric(runif(N)>0.5)
  
  obs_response <- ifelse(treatment==1,
                         as.numeric((P[2,1] * (1 - theta) + P[2,2]*theta)>runif(N)),
                         as.numeric((theta - bias)>runif(N)))
  
  out_data <- tibble(y = obs_response,
                     treatment=treatment)
  
  # define custom family as being upper and lower bounded for bias
  
  family_sens_reg <- custom_family("sens_reg",
                                 dpars=c("mu","bias"),
                                 links=c("logit","logit"),
                                 type="int",
                                 lb=c(NA,0),
                                 ub=c(NA,1),
                                 vars=c("P","vint1[n]"))
  
  # estimate model with brms

  est_mod_brms <- brm(y | vint(treatment) ~ 1,
                      family=family_sens_reg,
                      stanvars=all_stanvars,
                      data=out_data,
                      #prior=prior(normal(0,10),class="bias"),
                      prior=prior(beta(1,1),class="bias") +
                        prior(normal(-0.84,10),class="Intercept"),
                      chains=1,
                      cores=1,
                      iter=1000,
            backend = "cmdstanr")
  
  get_est_brms <- posterior::summarise_draws(as_draws(est_mod_brms, "b_Intercept"),
                                             estimate=~median(plogis(.x)),
                                             high=~quantile(plogis(.x),.95),
                                             low=~quantile(plogis(.x),.05)) %>% 
    mutate(param="mu") 
  
  # calculate bias transformed
  
  bias_est <- as_draws_matrix(est_mod_brms, "bias")
  
  mu_trans <- posterior_epred(est_mod_brms) %>% apply(1,median)

  bias_trans_est <- tibble(`95%`=quantile(bias_est * mu_trans,.95),
                         estimate=median(bias_est * mu_trans),
                         `5%`=quantile(bias_est * mu_trans,.05)) %>% 
    mutate(param="bias")
  
  get_est_brms <- bind_rows(get_est_brms, bias_trans_est) %>% 
    mutate(model="Bayes BRMS")
  
  # compare with RRreg
  
  est_freq <- RRuni(response=obs_response[treatment==1],
                  p=.25,
                  model="Crosswise",MLest = TRUE)
  
  this_sum <- summary(est_freq)
  
  # now do cheap Bayes
  
  lambda <- mean(obs_response[treatment==1])
  
  sim_data2 <- list(success_cm=(sum(treatment))*((lambda - .75)/(-.5)),
                     fail_cm= (sum(treatment))*(1-((lambda - .75)/-.5)),
                     success_d=sum(obs_response[treatment==0]),
                     fail_d=sum(1 - obs_response[treatment==0]))

  fit_simple <- sen_mod$sample(chains=1, cores=1, data=sim_data2,refresh=0)
  
  simple_sum <- fit_simple$summary()
  
  mu <- try(tibble(estimates=c(simple_sum$median[simple_sum$variable=="pi_hat"],
                     this_sum$coefficients[1,1],
                     get_est_brms$estimate[get_est_brms$param=="mu"]),
         q5=c(simple_sum$q5[simple_sum$variable=="pi_hat"],
                     this_sum$coefficients[1,1] - 1.96*this_sum$coefficients[1,2],
              get_est_brms$`5%`[get_est_brms$param=="mu"]),
         q95=c(simple_sum$q95[simple_sum$variable=="pi_hat"],
                     this_sum$coefficients[1,1] + 1.96*this_sum$coefficients[1,2],
               get_est_brms$`95%`[get_est_brms$param=="mu"]),
         models=c("RRreg","Bayes Simple","Bayes BRMS"),
         param="mu") %>% 
    mutate(cov=theta > q5 & theta < q95,
           sim=s))
  
  bias_est <- try(filter(get_est_brms, param=="bias") %>% 
    select(estimates="estimate",
           q5="5%",
           q95="95%") %>% 
    mutate(cov=bias > q5 & bias < q95,
           sim=s,
           param="bias",
           models="brms"))
  
  try(bind_rows(mu, bias_est))
   
  
}) %>% bind_rows
  
  saveRDS(over_sims, "data/over_sims.rds")
  
} else {
  
  over_sims <- readRDS("data/over_sims.rds")
  
}



over_sims %>% group_by(models,param) %>% 
  summarize(mean_cov=mean(cov)) %>% 
  knitr::kable(caption="Coverage")
Table 1: Coverage
models param mean_cov
Bayes BRMS mu 0.906
Bayes Simple mu 0.956
RRreg mu 0.568
brms bias 0.908
over_sims %>% group_by(models,param) %>% 
  summarize(spread_CIs=sum(q95 - q5)) %>% 
  knitr::kable(caption="Total Variance")
Table 1: Total Variance
models param spread_CIs
Bayes BRMS mu 71.44307
Bayes Simple mu 96.11747
RRreg mu 37.30712
brms bias 73.93706

Now let’s check bias.

#RMSE

over_sims %>% group_by(models,param) %>% 
  summarize(mean_rmse=switch(unique(param),
                             mu=sqrt(mean((estimates - theta)^2)),
                             bias=sqrt(mean((estimates - bias)^2)))) %>% 
  knitr::kable(caption = "RMSE")
Table 2: RMSE
models param mean_rmse
Bayes BRMS mu 0.0419459
Bayes Simple mu 0.0485183
RRreg mu 0.0487269
brms bias 0.0459833
over_sims %>% group_by(models,param) %>% 
  summarize(abs_bias=switch(unique(param),
                             mu=mean(abs(estimates - theta)),
                             bias=mean(abs(estimates - bias)))) %>% 
  knitr::kable(caption = "Mean Absolute Bias")
Table 2: Mean Absolute Bias
models param abs_bias
Bayes BRMS mu 0.0341473
Bayes Simple mu 0.0388322
RRreg mu 0.0390247
brms bias 0.0377768

References

Blair, Graeme, Kosuke Imai, and Yang-Yang Zhou. 2015. “Design and Analysis of the Randomized Response Technique.” Journal of the American Statistical Association 110 (511): 1304–19. https://doi.org/10.1080/01621459.2015.1050028.
Gingerich, Daniel W., Virginia Oliveros, Ana Corbacho, and Mauricio Ruiz-Vega. 2016. “When to Protect? Using the Crosswise Model to Integrate Protected and Direct Responses in Surveys of Sensitive Behavior.” Political Analysis 24 (2): 132–56. https://doi.org/10.1093/pan/mpv034.
Hout, Ardo van den, and Peter G. M. van der Heijden. 2002. “Randomized Response, Statistical Disclosure Control and Misclassification: A Review.” International Statistical Review / Revue Internationale de Statistique 70 (2): 269–88. https://doi.org/10.2307/1403910.
To leave a comment for the author, please follow the link and comment on their blog: R on Robert Kubinec.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Exit mobile version