# Duplicating stata survey design using svrepdesign from survey package in R

I am trying to use the R package 'survey' to produce some weighted descriptors of a dataset (Fragile Families) for a meta-analysis. According to the documentation, the weights are provided by the variable k4natwt, and the replication weights are provided by a series of variables: 'k4natwt_rep1 - k4natwt_rep26'.

The documentation further specifies the correct svyset command in Stata would be:

svyset [pweight=BASICWEIGHT], jkrw(REPLICATES, multiplier(1)) vce(jack) mse


I'm trying to replicate this survey design in R, but I think I am doing it incorrectly. I have a loose knowledge of the relevant statistics, and I suspect I am incorrectly translating the multiplier(1) argument in Stata in to the rscales argument in survey (for that matter, I may be incorrectly assuming that they are getting at the same statistical construct). Here's what I have in R so far (here I'm just trying to get a weighted mean of the cm4b.age variable).

> y5design<-svrepdesign(repweights="k4natwt.rep[1-9]+",type="JKn",weights=~k4natwt,data=y5weighted,combined.weights=TRUE,rscales=1)
> y5design
Call: svrepdesign.default(repweights = "k4natwt.rep[1-9]+", type = "JKn",
weights = ~k4natwt, data = y5weighted, combined.weights = TRUE,
rscales = 1)
Stratified cluster jackknife (JKn) with 26 replicates.
> svymean(y5weighted$cm4b.age,y5design,na.rm=TRUE) Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extent  So here are my questions: 1. Is this the most correct svrepdesign statement for the design? 2. What am doing wrong such that I get the dimnames error? ## 3 Answers Yes, scale=1, rscales=1 is correct. The general variance formula (from survey::svrVar) is  scale * sum( rscales*(replicates - centering) )  where replicates are the replicate values and centering is the mean of the replicates for mse=FALSE or the overall point estimate for mse=TRUE. The Stata syntax multiplier(1) is equivalent to scale=1, rscales=1. The Stata default is the same as the R default, with rscales equal to $$n_h/(n_h-1)$$ where $$n_h$$ is the number of replicates (and of PSUs) for stratum $$h$$. To specify non-default, non-constant rscales in Stata, you put a set of space-separated numbers in the multiplier() option. There is arguably a bug in svrepdesign in that you shouldn't have to specify scale=1 if you specify type="JKn" and rscales=. As it happens, if you don't specify the scale argument it stays as NULL and that's what causes the error you're seeing. replace y5weighted$cm4b.age with ~cm4b.age

• Hmm that does look correct but I still get the same error. May 22 '19 at 22:25
• It does resolve the error if I set both the scale and rscales arguments to 1. How can I tell if this is correct/what it ought to be? Is that suggested by the Stata syntax multiplier(1)? May 28 '19 at 17:04

I've discovered that setting both scale and rscales resolves the error, so I just need to set them correctly. I found a helpful forum post about what is meant by the svyset syntax. This suggests:

• In Stata, jkrweight will use 'jacknife n' by default, so that's consistent with the R code.
• vce(jack) just indicates that jacknife should be used, which is also consistent.
• 'mse' indicates that svyset uses the mean squared error formula for the variance calculation. We can specify this in svrepdesign using mse=TRUE (according to the svrepdesign documentation).
• According to the srvset documentation, multiplier indicates "the value of a jackknife multiplier to be added as a characteristic of the jackknife replicate-weight variables." This does appear to be consistent with the scale and rscales arguments to svrepdesign, save that svrepdesign allows both a global and replicate specific multiplier. Since the Fragile Families data specifies a value of 1, it seems to me that both scale and rscales sholud be 1. This is consistent with most of the examples I've found and produces sensible output, so that is likely the solution.