Using "Medsurvey" Package for Mediation analysis of complex survey data

[prest]

Hello here. I am new to R but I want to do mediation analysis on complex survey data using med.fit.brr in the medsurvey package. I did some reading about this method. The method requires use of Balanced Replication weights. In my data, I have the normal survey design elements (psu, strata and weight variables).

I tried to use this to create balanced replication weights;

meddesign<-svydesign(data=med,
                     id=~psu, strata=~kstrata,
                     weights=~pooled_wt2, nest=TRUE)
med.brr<- 
  as.svrepdesign( 
   meddesign , 
    type = "BRR" , 
    mse = TRUE ,
    compress = FALSE)

I got this error:-

# Error in brrweights(design$strata[, 1], design$cluster[, 1], ..., fay.rho = fay.rho,  : 
  # Can't split with odd numbers of PSUs in a stratum

Qn1. How can I address this?

To try to solve this, I used the Jackknife replication for stratified samples;

med.Jkn<- 
  as.svrepdesign( 
   meddesign , 
    type = "JKn" , 
    mse = TRUE ,
    compress = FALSE)

And weights were successfully created and named "repweights".
I then tried to use these weights to run mediation analysis using the med.fit.BRR;

# Specification of the model
model2<- '# outcome 
HE_glu~ u0*1 + c*diinormal + b*HE_hsCRP + age+ sex+ N_EN + DM_fh + HE_HP
# mediator 
HE_hsCRP~ u1*1 + a*diinormal +age +sex+ N_EN+HE_BMI+HE_HP+DM_fh  
# indirect effect (a*b) 
ab := a*b
# total effect 
total := c + (a*b) '

# Mediation analysis
Fit.brr<-med.fit.BRR(model=model2, data=med, 
                     mwgtname= pooled_wt2, 
                     repwgtnames= repweights, 
                     fayfactor=0.5, parallel='parallel', ncore=2)

I got this error

# Error in med.fit.BRR(model = model2, data = med, mwgtname = pooled_wt2, : value of HE_gluis not numerical

I am sorry it is a long procedure but I am kindly requesting for assistance.

[nirgrahamuk]

Seems, your key issue is HE_glu, the function med.fit.BRR wants it to be a numeric type.
Is it perhaps a factor ? would it have a natural translation to numeric that would work if you called as.numeric() on it ?

[StatSteph]

I agree that is the error (HE_gluis is not a numeric variable) but this function expects BRR weights and you have jacknife weights so it still won't work for the variances.

[prest]

Thank you so much.

The variable HE_glu is plasma glucose in mg/dL and it is therefore a continuous variable.

Do I still need to try calling it to numeric?

[prest]

Hello. Thank you for this insight. I suspected this, but I wanted to see how this works.

Is there a way I could create BRR weights with an odd number of psu in a stratum?

[nirgrahamuk]

Please post the result of running

 dput(head(med$HE_glu))

[prest]

@nirgrahamuk

dput(head(med$HE_glu))

Hello,
Below are the results.

dput(head(med.survey$HE_glu))
c(85, 86, 98, 91, 83, 79)

Also, kindly check the way I wrote the names of the replication weights in this code:
out.model<-med.fit.BRR(model=model.specify, data=med.survey,
mwgtname= pooled_wt2,
repwgtnames=(RepWt_1-RepWt_575),fayfactor=F,
parallel= NULL)

[nirgrahamuk]

wow, thats suprising about HE_glu, as it does seem numeric.
how big is your med data, too big to share?
dim(myd) will till row/column dimensions
and object.size(myd) will tell the byte size.

Is the data you pulled the same as the data we were discussing, because in your first shared codes the data was med, but now you are showing med.survey.
Are they the same, or different in some way ?

[prest]

dim(med.survey)
[1] 9058 602
object.size(med.survey)
44165992 bytes

They are somehow different. I exported it to sas and created Jacknife replication weights from there since I am somehow used to sas. I wanted to have one data set with replication weights and the other variables. They are 575 replication weights.

[prest]

Please find data.

med_rep.sas7bdat