Test Code

[Steve Samuels]

The link to the SEELS study that you provide does not work for me. This study had three waves. For a longitudinal weight, you've chosen the student weight from wave 3. In the absence of other information (that may be available to you), I would have chosen the average of weights for the three waves.

Below, I show code for doing BRR analysis "by hand". The idea is to create a different MEGLM command for each replicate and combine the results at the end.

You start by finding a good model, Do this with meglm that ignores BRR and uses the base weight wt_sa3. So, you will check for interactions and possible non-linear terms. Don't worry, you can check your model with the BRR analysis. Your current model is not correct. If you study the survey section of the meglm manual entry,you will see that you need a weight for every model level. The weight in the [pweight = ] option is for the level one units, in your case, waves. As all waves are included, you can set the first level weight to 1. The student weight should go in the student level section of the command.

Unfortunately, me models frequently fail for some replicates, and these failures can destroy the "balance" in "BRR" ("balanced repeated replication").

You currently use an if condition to define the subpopulation (treatment !=0). Doing can lead to incorrect standard errors. Create your own "subpop" option by setting the wave weight to 1 only for people in the subpopulation.

Here is my guess at code for your problem.

matrix drop _all  //house cleaning
gen wavwt = 1* (_treat !=0)   // change for other subpopulations or make equal to 1 for whole population
/* Step 1 Run a "good model"  */
meglm sa3lw_w age retain_k  [pw =wavwt], || StudentID:, cov(un) pweight(wt_sa3)
matrix b0 = e(b) //save the coefficients

/* To use a fovalues loop,r emove 0 from the  name of first nine weights. 
forvalues j = 1/9 {
    rename wt_sa3_0`i' wt_sa3`i'
}
/*  loop through the replicate weights, saving the results in a matrix. */
fovalues i = 1/32 {
    qui cap meglm sa3lw_w age retain_k  ///
   [pw = wavwt*(_treat!=0)], ///
   || StudentID:, cov(un) pweight(wt_sa3_0`i')
   if _rc  {
    di "replicate `i' failed"
    }
   matrix betas = nullmat(betas) / e(b)
   scalar rgood = rowsof(betas)  //count remaining replicates
 
/* Define a program to calculate the standard and root mean square errors of the betas matrix */
cap program drop _all
program define combine_reps, eclass
  mata: bb0 = st_matrix("b0")
  mata: bb = st_matrix("betas")
  mata: r = rows(bb)
  mata: v1 = (r-1)*variance(bb)/r
  mata: dif = bb:- st_matrix("b0")
  mata: v2 = cross(dif,dif)/r
  mata: st_matrix("v1",v1)
  mata: st_matrix("v2",v2)
  /*Display Standard errors */
  ereturn repost b = b0  V = v1
  di as text" Standard Errors based on "  as result rgood  " replicates"
  ereturn display
  mata: st_matrix("b0",bb0)  //b0 is lost
  ereturn repost b = b0  V = v2
  di as text " Root Mean Square Errors  "
  ereturn display
end
combine_reps  // run the program

Useful References:
variance estimation section of the Stata Survey Manual

Kolenikov, S. (2010). Resampling variance estimation for complex survey data. Stata Journal, 10(2), 165-199.
http://www.stata-journal.com/sjpdf.h...iclenum=st0187