DeclareDesign Community

Simulating Contrasts for Power Calculations

Let’s say I have the following design in which all estimates are informed by prior pilot data. Major thank you to @nfultz for holding my hand through the whole process here: Power calculations in non-standard design

ate <- 9.16 # size of coefficient on outgroup from from lm_robust(feelings_change ~ outgroup_pairing, data = data, team = team_id)

# Mean and Standard Deviation for the "pre" feelings measure for control condition
mean_pre_control = 19.6 
sd_pre_control = 22.5

# Mean and Standard Deviation for the "post" feelings measure for control condition
mean_post_control =  20.3 
sd_post_control = 23.0

# Mean and Standard Deviation for the "pre" feelings measure for treatment condition
mean_pre_treatment = 13.9
sd_pre_treatment = 18.6

# Mean and Standard Deviation for the "post" feelings measure for treatment condition
mean_post_treatment = 23.8
sd_post_treatment = 24.8

pop <- declare_population(N=600, 
                          party=rep(c('D','R'), each=N/2), 
                          
                          # generate draws from a truncated normal since these are "feelings thermometers"
                          # bounded between 0 and 100  
                          feelings_toward_outgroup_pre = msm::rtnorm(N, 
                                                                     mean = mean_pre_control, 
                                                                     sd = sd_pre_control,
                                                                     lower = 0, upper = 100),
                          feelings_toward_outgroup_post = msm::rtnorm(N, 
                                                                      mean = mean_post_control, 
                                                                      sd = sd_post_control,
                                                                      lower = 0, upper = 100)
)

# Create an "assn_team" function that generates random pairings of DD, RR, RD, and DR
assn_team <- declare_assignment(handler=function(data) {
  N <- nrow(data)
  D_teams <- c(paste("D", rep(1:(N/8), times=2)), paste("X", 1:(N/4)))
  R_teams <- c(paste("R", rep(1:(N/8), times=2)), paste("X", 1:(N/4)))
  
  data$team <- NA
  data$team[data$party == 'D'] <- D_teams
  data$team[data$party == 'R'] <- R_teams
  
  data$Z <- +grepl("X", data$team)
  
  # Variable that gets the party of one's partner (sanity check)
  data$partner_party =  with(data, ave(party, team, FUN = rev))
  
  data
})

df <- assn_team(pop())

potential_outcomes <- declare_potential_outcomes(Y_Z_0 = feelings_toward_outgroup_post - feelings_toward_outgroup_pre, 
                                                Y_Z_1 = (feelings_toward_outgroup_post - feelings_toward_outgroup_pre) + ate)

df <- potential_outcomes(df)

estimand <- declare_estimand(ATE = mean(Y_Z_1 - Y_Z_0))
reveal <- declare_reveal(Y, Z)
estimator <- declare_estimator(Y ~  party * Z, 
                               estimand = estimand, 
                               model = lm_robust, 
                               clusters = team,
                               term = "Z")
design <- pop + assn_team + potential_outcomes+ estimand  + reveal + estimator
diagnosis <- diagnose_design(design)

This results in the following output:

   Design Label Estimand Label Estimator Label Term N Sims   Bias   RMSE  Power Coverage Mean Estimate SD Estimate Mean Se Type S Rate Mean Estimand
         design            ATE       estimator    Z    500   0.12   2.90   0.90     0.94          9.28        2.90    2.86        0.00          9.16
                                                           (0.12) (0.09) (0.01)   (0.01)        (0.12)      (0.09)  (0.01)      (0.00)        (0.00)

In my actual study, I conduct contrasts like so:

 model = lm_robust(feelings_change ~ politics * outgroup_pairing, 
                 cluster = team_id,
                 se = "stata",
                 data = data)
 rg = qdrg(object = model, data = data)
  outgroup_means = emmeans(rg, ~ outgroup_pairing | politics)
  contrast_text = contrast(outgroup_means, method = "pairwise") 
politics = Democrat:
 contrast estimate   SE  df t.ratio p.value
 0 - 1      -11.66 3.33 199 -3.502  0.0006 

politics = Republican:
 contrast estimate   SE  df t.ratio p.value
 0 - 1       -6.66 3.27 199 -2.040  0.0427 

I want to make sure that I’m sufficiently powered to be able to conduct these types of contrasts, but cannot decide the optimal way to do so within DeclareDesign, because contrasts are conducted with respect to the model itself. One idea is to re-parametrize so that there are four “treatment” conditions (DD, DR, RD, and RR), but I’m not sure how to make explicit the comparisons between (DD v DR) and (RR vs RD).

Thanks for any assistance in advance.

One random question: I noticed that when I don’t cluster at the team-level, my power actually goes down. This is a bit strange, since I assumed clustering would nearly always inflate standard errors. Is this strange behavior?

For linear contrasts, see ?estimatr::lh_robust - for nonlinear contrasts, I had written a section in the postestimation vignette.

I think for your case, lh_robust would suffice because you had a balanced design and match the marginal, but double check it against emmeans to be sure. If there’s a big difference, writing a wrapper function for emmeans shouldn’t be too difficult.

This is my intuition as well, so it is a bit odd, but I think it’s possible for the uncorrected variance to be biased in either direction.

Also if you were to use mixed effects models instead of robust standard errors to deal with the clustering, you would expect power to increase because the variance from clustering gets partitioned out.

Also https://stats.stackexchange.com/questions/120507/cluster-robust-standard-errors-are-smaller-than-unclustered-ones-in-fgls-with-cl

EDIT:
which links to this, which seems relevant:

https://www.statalist.org/forums/forum/general-stata-discussion/general/329139-cluster-robust-standard-errors-are-smaller-than-unclustered-ones-in-fgls-with-cluster-fixed-effects

Thanks so much @nfultz. I think I’m confused as to how I’m supposed to use lh_robust (I’ve had problems with this in the past, which is why I figured out how to get emmeans to work with lm_robust.?

  1. For the specific contrast noted above, in which I’m comparing means of outgroup (a 0 or 1 variable) at different levels of party, how would I do that with a single linear hypothesis specification, as required in lh_robust

  2. I’m having a more fundamental confusion about where in the declare-design order of operations a contrast would go (either from an emmeans wrapper or lh_robust). Does it go in the declare_estimator stage?

Thanks for the links about clustered standard errors. It looks like clustered standard errors can increase or decrease standard errors depending on the direction of the ICC.


For example, using emmeans and follow your tutorial, I’ve tried to hack together something like this:

contrast_handler <- function(data) {
  model = lm_robust(Y~ as.factor(Z) * party, data = data,
            clusters = team)
  rg = qdrg(object = model, data = data)
  means = emmeans(rg, ~ Z | party)
  contrast_text = contrast(means, method = "pairwise")
  #print(contrast_text)
  tidy(contrast_text) %>% unite(party, level1:party)
}

design_og <- pop + assn_team + potential_outcomes + estimand  + reveal + contrast_estimator
diagnosis <- diagnose_design(design_og)

This seems to run okay, with the following results:

Research design diagnosis based on 1000 simulations. Diagnosand estimates with bootstrapped standard errors in parentheses (100 replicates).

 Design Label Estimand Label Estimator Label N Sims   Bias   RMSE  Power Coverage Mean Estimate SD Estimate Mean Se Type S Rate Mean Estimand
    design_og            ATE       estimator   1000  -0.14   0.14   0.71       NA         -0.07        0.03    0.03        1.00          0.07
                                                    (0.00) (0.00) (0.01)       NA        (0.00)      (0.00)  (0.00)      (0.00)        (0.00)

However, I don’t know what the estimated power of .71 applies to. There are two sets of contrasts… does the .71 power apply to the ability to detect both of the contrasts?

There’s some edge cases around custom estimators and integrating with estimands. Running your handler directly definitely shows the two results:

> draw_data(design) %>% contrast_handler()
# A tibble: 2 x 6
  party estimate std.error    df statistic   p.value
  <chr>    <dbl>     <dbl> <dbl>     <dbl>     <dbl>
1 0_1_D    -8.69      2.96   596     -2.94 0.00342  
2 0_1_R   -12.2       2.87   596     -4.23 0.0000267

^^^ Also note that you get 0 - 1 there.

Since you have two estimates, I would suggest having two different estimands, one for each party, and make them match the estimator.

If you provide them in order to a tidy_estimator step, they will be joined in order as well:

contrast_handler <- function(data) {
  model = lm_robust(Y~ as.factor(Z) * party, data = data,
                    clusters = team)
  rg = qdrg(object = model, data = data)
  means = emmeans(rg, ~ Z | party)
  contrast_text = contrast(means, method = "pairwise")
  #print(contrast_text)
  tidy(contrast_text) %>% unite(term, level1:party)
}

estimand_D <- declare_estimand(ATE_D=mean(Y_Z_0 - Y_Z_1), subset= party == "D")
estimand_R <- declare_estimand(ATE_R=mean(Y_Z_0 - Y_Z_1), subset= party == "R")

contrast_estimator <- declare_estimator(handler = tidy_estimator(contrast_handler), estimand=c("ATE_D", "ATE_R"))

design_og <- pop + assn_team + potential_outcomes + estimand_D + estimand_R  + reveal + contrast_estimator


diagnosis <- diagnose_design(design_og)

Which yields:

 Design Label Estimand Label Estimator Label  Term N Sims   Bias   RMSE  Power Coverage Mean Estimate SD Estimate Mean Se Type S Rate Mean Estimand
    design_og          ATE_D       estimator 0_1_D    500   0.04   2.82   0.88       NA         -9.12        2.82    2.84        0.00         -9.16
                                                          (0.12) (0.09) (0.01)       NA        (0.12)      (0.09)  (0.01)      (0.00)        (0.00)
    design_og          ATE_R       estimator 0_1_R    500  -0.02   2.85   0.89       NA         -9.18        2.86    2.85        0.00         -9.16
                                                          (0.11) (0.09) (0.01)       NA        (0.11)      (0.09)  (0.01)      (0.00)        (0.00)

If you go the lh_robust route, you use the regular estimator handler, and specify model=lh_robust and linear_hypothesis=''

estimator <- declare_estimator(Y ~  party * Z, 
                               estimand = estimand, 
                               model = lh_robust, 
                               clusters = team,
                               linear_hypothesis=c("Z=0", "partyR:Z + Z = 0"),
                               term=TRUE)

which picks up the LH terms:

> draw_data(design) %>% estimator
  estimator_label             term      estimate   std.error      statistic         p.value     conf.low    conf.high          df outcome estimand_label
1       estimator      (Intercept) -0.2147655956 2.004175472 -0.10715907791 0.9149527742938 -4.208172396  3.778641205  74.0000000       Y            ATE
2       estimator           partyR -0.2693402441 2.935806497 -0.09174318689 0.9270261440701 -6.070853500  5.532173011 148.0000000       Y            ATE
3       estimator                Z 10.1903799058 2.936153236  3.47065670186 0.0006373786853  4.400194043 15.980565769 197.7757848       Y            ATE
4       estimator         partyR:Z -0.2064312013 3.989831011 -0.05173933449 0.9587711760552 -8.058353265  7.645490862 296.9966330       Y            ATE
5       estimator              Z=0 10.1903799058 2.936153236  3.47065670186 0.0005567275216  4.423915113 15.956844699 596.0000000       Y            ATE
6       estimator partyR:Z + Z = 0  9.9839487046 2.906115799  3.43549582906 0.0006326037521  4.276476005 15.691421405 596.0000000       Y            ATE

Thanks so much for this! This is immensely helpful.