You can try using the following option which does a conservative adjustment of adjusting the lonely-PSU (a stratum with only one PSU) around the grand mean. I'm not sure why this wouldn't occur in the interactive session but does occur in markdown unless you applied this option and forgot about it.
options(survey.lonely.psu="adjust")
Now, that I look at your code, I see you try to put a survey.lonely.psu option when you create the design object but that's not where it goes. There's some other odd things in your code. I've redone it below correctly using both survey package and the srvyr package. I'm just doing analysis for Italy to make it run a bit slower as an example.
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.3 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(survey)
## Loading required package: grid
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Loading required package: survival
##
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
##
## dotchart
library(here)
## here() starts at C:/Users/..
# ESS data here: https://www.dropbox.com/s/ceq02i2rkyaqrwk/ESS9e03_1.sav
# brief codebook https://www.dropbox.com/s/d9jefeuht8bu6yo/appendix%20general%20codebook.pdf
ess9_in <-haven::read_spss(here::here("ESS9e03_1.sav"))
ess9<-ess9_in %>%
mutate(
age=agea,
age = na_if(agea, 999),
agecat= case_when(
age %in% 15:25 ~ "1. 15 to 25",
age %in% 26:36 ~ "2. 26 to 36",
age %in% 37:47 ~ "3. 37 to 47",
age %in% 48:58 ~ "4. 48 to 58",
age %in% 59:69 ~ "5. 59 to 69",
age %in% 70:90 ~ "6. 70 to 90"),
agecat = fct_relevel(agecat, "1. 15 to 25","2. 26 to 36","3. 37 to 47","4. 48 to 58","5. 59 to 69","6. 70 to 90"),
trust_parliament = as.factor(trstprl),
tp_num=as.character(trust_parliament) %>% as.numeric(),
anweight=pweight*pspwght
) %>%
select(psu, stratum, anweight, agecat, trust_parliament, tp_num, cntry)%>%
filter(cntry=="IT")
options(survey.lonely.psu="adjust")
ess_design<-svydesign(weights=~anweight, ids=~psu, strata=~stratum, data = ess9, nest=TRUE)
svyby(~tp_num, ~ agecat, ess_design, na.rm=TRUE, svymean)
## agecat tp_num se
## 1. 15 to 25 1. 15 to 25 4.640104 0.1577266
## 2. 26 to 36 2. 26 to 36 4.279145 0.1589022
## 3. 37 to 47 3. 37 to 47 4.364255 0.1460727
## 4. 48 to 58 4. 48 to 58 4.413230 0.1417483
## 5. 59 to 69 5. 59 to 69 4.103305 0.1499386
## 6. 70 to 90 6. 70 to 90 3.778527 0.1466672
library(srvyr)
##
## Attaching package: 'srvyr'
## The following object is masked from 'package:stats':
##
## filter
ess_design_2 <- ess9 %>%
as_survey_design(ids=psu, strata=stratum, weights=anweight, nest=TRUE)
ess_design_2 %>%
group_by(agecat) %>%
summarise(
tp=survey_mean(tp_num, na.rm=TRUE)
)
## # A tibble: 7 x 3
## agecat tp tp_se
## <fct> <dbl> <dbl>
## 1 1. 15 to 25 4.64 0.158
## 2 2. 26 to 36 4.28 0.159
## 3 3. 37 to 47 4.36 0.146
## 4 4. 48 to 58 4.41 0.142
## 5 5. 59 to 69 4.10 0.150
## 6 6. 70 to 90 3.78 0.147
## 7 <NA> 3.85 0.374