1. Load in the needed functions and data

First, load in the mixPV function. Load in the pisa2012 dataset as well.

library(modelsummary)
source("https://raw.githubusercontent.com/flh3/pubdata/main/mixPV/mixPVv2.R")
data(pisa2012, package = 'MLMusingR') 

By default, the WeMix package uses conditional weights = FALSE or cWeights = FALSE (no need to specify). This makes it easier so users don’t have to compute the conditional weights (conditional on the school being selected) and can just use the total weight supplied in many datasets. You have to make sure which type of weight you are using. NOTE: the mc = TRUE option is used to shorten computing time (turning on the use of multiple cores). The following just uses the raw weights at both levels:

m1a <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               sc14q02 + st04q01 + escs + (1|schoolid), 
             weights = c('w_fstuwt', 'w_fschwt'), mc = TRUE,
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

However, if you use software such as Mplus or SAS, you need to specify the conditional weights. To compute the conditional weights (if you are specifying weights at both levels), it is the total weight divided by the school weight.

pisa2012$cweight <- pisa2012$w_fstuwt / pisa2012$w_fschwt

We can then use this when we fit the model and we include cWeights = TRUE. We can compare the output in a while.

m1b <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               sc14q02 + st04q01 + escs + (1|schoolid), 
             weights = c('cweight', 'w_fschwt'), mc = TRUE,
             cWeight = TRUE,
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

I will include the example if we just use the conditional weight (as you would with other software packages) but DO NOT turn on the cWeights option.

m1c <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               sc14q02 + st04q01 + escs + (1|schoolid), 
             weights = c('cweight', 'w_fschwt'), mc = TRUE,
             #cWeight = TRUE,  ## for illustrative purposes only
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

We can compare results side-by-side.

modelsummary(list("default" = m1a, "cweights" = m1b, "wrong" = m1c),
  stars = TRUE, gof_map = NA)

The first two models have identical results. The third model (which is incorrectly specified) results differ and the level-2 variance here is zero.

2. What about scaling level-1 weights?

Asparouhov (2005) suggested scaling the weights at level 1. The following formulas (and explanations) can be seen at the Mplus website at https://statmodel.com/download/Scaling3.pdf.

The unscaled weight is \(w^*_{ij}=w_{ij}\).

What Mplus calls the cluster weight has scaled weights that add up to the cluster sample size. Carle (2009) refers to this as Method A. It is computed using:

\(w^*_{ij}=w_{ij}\frac{n_j}{\Sigma_i w_{ij}}\)

Here, \(n_j\) is the number of sample units in cluster \(j\). This is the default scaling method in Mplus and has generally been shown to perform well (see Asparouhov, 2006).

Another scaling approach is referred to as the ecluster method which is the effective cluster sample size method. This is referred to as Method B by Carle (2009) (though the formula in his appendix seems to differ slightly):

This is similar to the cluster method except that the numerator is \(n^*_j\) instead and is:

\(w^*_{ij}=w_{ij}\frac{n^*_j}{\Sigma_i w_{ij}}\).

The numerator is computed using:

\(n^*_{j}=\frac{(\Sigma_i w_{ij})^2}{\Sigma_i w^2_{ij}}\).

To compute the scaled weights, we can use the wscale function (note the options). The cluster variable, dataset, conditional cluster weight, and type of scaling need to be specified:

pisa2012$cl_weight <- wscale(cluster = 'schoolid', data = pisa2012, 
                             wt = 'cweight', type = 'cluster')
pisa2012$es_weight <- wscale(cluster = 'schoolid', data = pisa2012, 
                             wt = 'cweight', type = 'ecluster')

After I wrote the function above, I found a function that does the scaling. The results are the same:

xx <- datawizard::rescale_weights(pisa2012, "schoolid", "w_fstuwt")

We can compare results to those of Mplus. This is using cluster weights:

ex <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
              1 + (1|schoolid), 
            weights = c('cl_weight', 'w_fschwt'), mc = TRUE,
            cWeight = TRUE,
            data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

summary(ex)

Results of multilevel analyses with 5 plausible values.
Number of observations: 3136 

Estimates for random effects: 
                     estimate std.error statistic     df Pr(>t)    
schoolid.(Intercept)  1825.03    436.26      4.18 2719.0 <2e-16 ***
Residual              5726.58    233.29     24.55   90.6 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates for fixed effects: 
            estimate std.error statistic   df Pr(>t)    
(Intercept)   473.39      7.03     67.33 1777 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparing to Mplus output shows similar results (SEs might differ slightly):

                Estimate       S.E.  Est./S.E.    P-Value

Within Level

Variances
MATH            5726.580    234.000     24.473      0.000

Between Level

Means
MATH             473.390      7.387     64.088      0.000

Variances
MATH            1825.021    464.365      3.930      0.000

This is using the ecluster option:

ex2 <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               1 + (1|schoolid), 
             weights = c('es_weight', 'w_fschwt'), mc = TRUE,
             cWeight = TRUE,
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

summary(ex2)

Results of multilevel analyses with 5 plausible values.
Number of observations: 3136 

Estimates for random effects: 
                     estimate std.error statistic     df Pr(>t)    
schoolid.(Intercept)  1821.40    434.62      4.19 2714.0 <2e-16 ***
Residual              5727.56    233.36     24.54   90.2 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates for fixed effects: 
            estimate std.error statistic   df Pr(>t)    
(Intercept)   473.40      7.02     67.39 1776 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                Estimate       S.E.  Est./S.E.    P-Value

Within Level

Variances
MATH            5727.560    234.056     24.471      0.000

Between Level

Means
MATH             473.403      7.379     64.156      0.000

Variances
MATH            1821.389    462.665      3.937      0.000

This is unscaled:

ex3 <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               1 + (1|schoolid), 
             weights = c('cweight', 'w_fschwt'), mc = TRUE,
             cWeight = TRUE,
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

summary(ex3)

Results of multilevel analyses with 5 plausible values.
Number of observations: 3136 

Estimates for random effects: 
                     estimate std.error statistic   df Pr(>t)    
schoolid.(Intercept)  2126.22    503.39      4.22 3099 <2e-16 ***
Residual              5895.85    175.05     33.68  876 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates for fixed effects: 
            estimate std.error statistic   df Pr(>t)    
(Intercept)   472.35      7.33     64.44 1470 <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                 Estimate       S.E.  Est./S.E.    P-Value
Within Level

Variances
MATH            5895.849    174.758     33.737      0.000

Between Level

Means
MATH             472.349      7.750     60.950      0.000

Variances
MATH            2126.477    535.991      3.967      0.000

Results differ based on the scaling option and results are all similar to those generated by Mplus.

3. What about just using level-2 weights?

A simpler option may just be to use level-2 weights (as suggested by Mang et al., 2021 as well as our own simulation results [under review]). To do so, when using R, we need to need to set the weight at level 1 to 1 and use conditional weights.

pisa2012$one <- 1 #level-1 weight
m1d <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               sc14q02 + st04q01 + escs + (1|schoolid), 
             weights = c('one', 'w_fschwt'), mc = T,
             cWeight = TRUE,
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

For comparison, I will show results also using the scaled weights:

m1e <- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
               sc14q02 + st04q01 + escs + (1|schoolid), 
             weights = c('cl_weight', 'w_fschwt'), mc = TRUE,
             cWeight = TRUE,
             data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

m1f<- mixPV(pv1math + pv2math + pv3math + pv4math + pv5math ~ 
              sc14q02 + st04q01 + escs + (1|schoolid), 
            weights = c('es_weight', 'w_fschwt'), mc = TRUE,
            cWeight = TRUE,
            data = pisa2012)

Attempting to use 15 cores. Progress will not be displayed. Please wait.

modelsummary(list("l2only" = m1d, 'cluster' = m1e, "ecluster" = m1f), 
             stars = TRUE, gof_map = NA)

Results are very similar.

A reason for this is that the level-1 scaled weights are all very close to 1 (or even equal to 1). If there is no variation within cluster, the weights will be 1.

If we examine the original conditional weights, there is quite some variability:

range(pisa2012$cweight) #the original conditional weight

[1]  1.228571 37.018006

However, compared to the scaled weights, they are all very close to 1. Which is why the results shown above are all very similar.

range(pisa2012$es_weight)

[1] 0.7635826 1.3325292

range(pisa2012$cl_weight)

[1] 0.7820465 1.3559951

mean(pisa2012$es_weight)

[1] 0.9962017

mean(pisa2012$cl_weight)

[1] 1

Takeaway: when running multilevel models with weights at different levels, make sure you are specifying the correct weights. If the residual variance at the second level is close to zero, you should probably revisit your syntax.