Inconsistent handling of lonely PSUs in R's `survey` package? There appears to be some inconsistent handling of lonely PSUs in R's survey package when calling svyglm on a subset of stratified data where the act of subsetting creates the lonely PSU. It appears that the sample size associated with a lonely PSU in a stratum does not change when using the subset.survey.design method, which means that svyglm runs without throwing an error indicating that the lonely PSU exists. Is it correct for svyglm to do so in that circumstance?
To further demonstrate the issue, I first create a dummy survey dataset.
library(survey)
library(tidyverse)

set.seed(9001)

# create random weights and rescale to add up to 1
wts = runif(1000, 1, 1000)
wts = wts/sum(wts)

# create dummy survey dataset for testing lonely PSU subsets
test_svy_data <- tibble(id = 1:1000, a = rnorm(1000), b = rnorm(1000), 
                        group = c(rep("a",999), "b"), wt = wts, 
                        strat = c(rep(1,333),rep(2,333),
                                  rep(3,332), c(4,4)))

This creates a dataset where subsetting to group == "a" creates a lonely PSU in stratum 4.
I then construct two survey.design objects using different ways of subsetting: 1) Filter the initial dataset so that the lonely PSU already exists, then create a survey.design object; and 2) Use the subset.survey.design method built into the survey package after creating a survey.design object, thereby producing the lonely PSU.
# create stratified survey from initial data with lonely PSU
test_svy1 <- svydesign(data = test_svy_data %>% filter(group == "a"),
                       ids = ~id, weights = ~wt, strata = ~strat)

# create stratified survey from initial data without lonely PSU, then subset
test_svy2 <- svydesign(data = test_svy_data,
                       ids = ~id, weights = ~wt, 
                       strata = ~strat) %>%
  subset(group == "a")

Next, I attempt to run svyglm using both of these test objects. The first test produces the expected lonely PSU error as per the survey package's documentation on lonely PSUs when using the default handling.
> try(test1 <- svyglm(a~b, test_svy1) %>% summary())
Error in onestrat(x[index, , drop = FALSE], clusters[index], nPSU[index][1],  : 
  Stratum (4) has only one PSU at stage 1

The second runs without any errors.
> try(test2 <- svyglm(a~b, test_svy2) %>% summary())
> # no error

The difference appears to arise from the fpc$sampsize component of the survey.design objects. In the first test, the sample size associated with each PSU is correct, per the following code and output. A sample size of one is associated with the lonely PSU.
> test_svy1$fpc$sampsize %>% table()
.
  1 332 333 
  1 332 666 

But after using subset.survey.design, the fpc$sampsize component still has a sample size of 2 associated with the lonely PSU.
> test_svy2$fpc$sampsize %>% table() 
.
  2 332 333 
  1 332 666
# appears to still have a sample size of 2 assigned to the lonely PSU

Assigning a value of 1 to the fpc$sampsize component of the survey.design object for the lonely PSU, then calling svglm using the new values restores the error:
> # try replacing fpc attribute with one for the lonely PSU
> test_svy3 <- test_svy2
> 
> # Assign corrected(?) fpc value
> test_svy3$fpc$sampsize[999,1] <- 1
> 
> # throws lonely PSU error
> try(test3 <- svyglm(a~b, test_svy3) %>% summary())
Error in onestrat(x[index, , drop = FALSE], clusters[index], nPSU[index][1],  : 
  Stratum (4) has only one PSU at stage 1

Is the package handling the lonely PSU created by a subset correctly in the first test? Or the second test? Both somehow?
[I have compared this handling to Stata's svyset approach. Stata warns about a lonely PSU regardless of when the subsetting occurs.]
Also for reference, here's the output from packageDescription("survey") run on my machine.
Package: survey
Title: Analysis of Complex Survey Samples
Description: Summary statistics, two-sample tests, rank tests, generalised linear models,
          cumulative link models, Cox models, loglinear models, and general maximum
          pseudolikelihood estimation for multistage stratified, cluster-sampled,
          unequally weighted survey samples. Variances by Taylor series linearisation or
          replicate weights. Post-stratification, calibration, and raking. Two-phase
          subsampling designs. Graphics. PPS sampling without replacement.
Version: 4.1-1
Author: Thomas Lumley
Maintainer: "Thomas Lumley" <t.lumley@auckland.ac.nz>
License: GPL-2 | GPL-3
Depends: R (>= 3.5.0), grid, methods, Matrix, survival
Imports: stats, graphics, splines, lattice, minqa, numDeriv, mitools (>= 2.4)
Suggests: foreign, MASS, KernSmooth, hexbin, RSQLite, quantreg, parallel, CompQuadForm,
          DBI, AER
URL: http://r-survey.r-forge.r-project.org/survey/
NeedsCompilation: no
Packaged: 2021-07-16 23:50:34 UTC; tlum005
Repository: CRAN
Date/Publication: 2021-07-19 08:40:04 UTC
Built: R 3.6.0; ; 2022-09-07 21:44:31 UTC; unix


There's a difference between lonely PSUs created by design (i.e., when a stratum only has one PSU sampled) and lonely PSUs created by subsetting to a domain (referred to as "domain lonely PSUs" in the terminology of the survey package).
Unbiased variance estimation for design-lonely PSUs is not generally possible for design-lonely PSUs, but for domain-lonely PSUs, variance estimation for totals is still unbiased. This is why the survey package quite sensibly doesn't throw errors for domain lonely PSUs but does throw errors for design-lonely PSUs.
If you're using the option(survey.lonely.psu = 'adjust') to handle design-lonely PSUs, you can use options(survey.adjust.domain.lonely = TRUE) to have this special handling used for domain-lonely PSUs as well. However, at the present moment, some caution is warranted around using option(survey.lonely.psu = 'adjust') as it has a bug when estimating population totals.. I think this will get fixed soon and I'll update this answer when it does.
The handling of lonely PSUs is documented in the survey package in a help file accessible in R by running help("survey.lonely.psu", package = "survey").
