I was comparing results that I generated in R for complex survey analysis using the survey package to results from SPSS using the complex samples analysis add-on. The sample size is large ~ N=5500

This is the R code:

svy <- svydesign(ids = ~cl1+houseID, strata = ~strata_study, weights = ~wt,  data=data_in,nest = TRUE, fpc = NULL)

svyby(~X, by = ~sadmood, design=svy, FUN=svymean)

With the output:

          X    group0    group1  se.group0  se.group1
No       No 0.7898876 0.2101124 0.03106912 0.03106912
Yes     Yes 0.7533348 0.2466652 0.04133818 0.04133818

In SPSS I navigated to the complex samples Analysis Preparation Wizard, selected sampling with replacement (WR) and the following:

Strata: strata_study

Clusters: cl1 (cluster 1), houseID (cluster 2)

Sample Weight: wt

The output as follows:

          X    group0    group1  se.group0  se.group1
No       No     79.0%     21.0%       1.4%       1.4%
Yes     Yes     75.3%     24.7%       4.0%       4.0%

While the estimates for % of group==0 and group==1 within X==NO and X==YES are similar in SPSS and R; the standard errors are different. Does anyone know why these differences are there?


  • 1
    $\begingroup$ It looks like you're using clusters in SPSS and strata in R (but I"m not v. familiar with complex samples). How many clusters / strata do you have. $\endgroup$ Mar 10, 2023 at 17:51
  • 1
    $\begingroup$ Thanks very much for the comment! This got me thinking about how R is specifying the clusters. I have 2 clusters (houseid which is nested within cl1). In the docmentation the smaller level should be presented second (it says id : "Formula or data frame specifying cluster ids from largest level to smallest level"), so I had cl1+houseid but now doubting whether this is right, as switching these gives me near identical results in R and SPSS! I have a feeling it should be houseid+cl1, unless its my SPSS specification thats incorrect $\endgroup$
    – s.stats
    Mar 10, 2023 at 18:34
  • 2
    $\begingroup$ ~cl1+houseid is correct for R. One way to check the SPSS is that dropping houseid from the specification will have little or no impact if the specification is correct, but could make a noticeable difference if it's wrong. [That is, 'with replacement' std errors only use the coarsest level of sampling] $\endgroup$ Mar 10, 2023 at 20:19
  • 1
    $\begingroup$ Is it possible that you specified both cluster variables at stage 1 of sampling? It appears from online documentation (ibm.com/docs/en/spss-statistics/…) that this is allowed and results in using all combinations of the variables to define PSUs. That would match with using ~houseid+cl1 in R, and would be wrong. $\endgroup$ Mar 10, 2023 at 20:31
  • $\begingroup$ Thanks @ThomasLumley for the comment! Yes both variables were specified at stage 1 in SPSS. Using any variation of cl1 and houseid , and when dropping cl1 in SPSS made no difference. The cluster size is very small ~1.07, so the effect of including houseid is negligible, there was no difference in the output for cl1+houseid or cl1 by itself in R. Similarly, the incorrect specification of houseid+cl1 in R gives the same output as houseid by itself; it either takes the first level only or appears this way possibly because the cluster size is so small. $\endgroup$
    – s.stats
    Mar 12, 2023 at 3:31

1 Answer 1


The R syntax is correct: ~cl1+houseid specifies that cl1 values identify sampling units at stage 1 (PSUs) and houseid values identify sampling units at stage 2. If instead you want to use combinations of two variables to identify PSUs, you need to use the interaction function to create a single variable with all combinations, eg, ~interaction(cl1, houseid). Specifying the formula backwards ~houseid+cl1 gives you houseid as the PSU (which, since they're nested, is the same as interaction(houseid, cl1).

In SPSS, you have the same options. At stage 1, specify just the PSU (cl1), and at stage 2 specify just the stage 2 sampling unit (houseid). Or, since stage 2 doesn't matter for 'with-replacement' standard errors, just specify stage 1. If you specify two cluster variables at stage 1 (as you did) you get all combinations of them as the PSUs, which is wrong.

There's extra potential for confusion in the social sciences because multilevel modellers call the finest partition of the data 'level 1' and survey samplers call the coarsest partition 'stage 1'.

Finally, I note that R will tell you how many clusters you have if you just print the survey design object. Here's one of the built-in examples, constructed the right way around and constructed backwards. The right way around, it describes itself correctly as a 2-stage sampling design and gives the numbers of clusters. The wrong way around, it describes itself as 'independent sampling' because we've told it that the PSUs come from what's really the second-stage identifier (snum), which identifies individual records in the dataset.

> dclus2
2 - level Cluster Sampling design
With (40, 126) clusters.
dclus2<-svydesign(id=~dnum+snum, fpc=~fpc1+fpc2, data=apiclus2)
> dclus2a<-svydesign(id=~snum+dnum, data=apiclus2,weight=~pw)
> dclus2a
Independent Sampling design (with replacement)
svydesign(id = ~snum + dnum, data = apiclus2, weight = ~pw)

I assume you can also get that sort of information in SPSS, but I don't speak SPSS.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.