Analysis

Data

Provided by {srvyrexploR} library. The first is America’s elections data.

Code

data('anes_2020', package = "srvyrexploR")
anes_2020 %>% select(-matches("^V\\d")) %>%# starts with V followed by a digit
glimpse()

Rows: 7,453
Columns: 21
$ CaseID                  <dbl> 200015, 200022, 200039, 200046, 200053, 200060…
$ InterviewMode           <fct> Web, Web, Web, Web, Web, Web, Web, Web, Web, W…
$ Weight                  <dbl> 1.0057375, 1.1634731, 0.7686811, 0.5210195, 0.…
$ VarUnit                 <fct> 2, 2, 1, 2, 1, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2…
$ Stratum                 <fct> 9, 26, 41, 29, 23, 37, 7, 37, 32, 41, 22, 7, 3…
$ CampaignInterest        <fct> Somewhat interested, Not much interested, Some…
$ EarlyVote2020           <fct> NA, NA, NA, NA, NA, NA, NA, NA, Yes, NA, NA, N…
$ VotedPres2016           <fct> Yes, Yes, Yes, Yes, Yes, No, Yes, No, Yes, Yes…
$ VotedPres2016_selection <fct> Trump, Other, Clinton, Clinton, Trump, NA, Oth…
$ PartyID                 <fct> Strong republican, Independent, Independent-de…
$ TrustGovernment         <fct> Never, Never, Some of the time, About half the…
$ TrustPeople             <fct> About half the time, Some of the time, Some of…
$ Age                     <dbl> 46, 37, 40, 41, 72, 71, 37, 45, 70, 43, 37, 55…
$ AgeGroup                <fct> 40-49, 30-39, 40-49, 40-49, 70 or older, 70 or…
$ Education               <fct> Bachelor's, Post HS, High school, Post HS, Gra…
$ RaceEth                 <fct> "Hispanic", "Asian, NH/PI", "White", "Asian, N…
$ Gender                  <fct> Male, Female, Female, Male, Male, Female, Fema…
$ Income                  <fct> "$175,000-249,999", "$70,000-74,999", "$100,00…
$ Income7                 <fct> $125k or more, $60k to < 80k, $100k to < 125k,…
$ VotedPres2020           <fct> NA, Yes, Yes, Yes, Yes, Yes, Yes, NA, Yes, Yes…
$ VotedPres2020_selection <fct> NA, Other, Biden, Biden, Trump, Biden, Trump, …

Code

anes_2020 %>% select(-matches("^V\\d")) %>% skimr::skim()

Data summary
Name	Piped data
Number of rows	7453
Number of columns	21
_______________________
Column type frequency:
factor	18
numeric	3
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
InterviewMode	0	1.00	FALSE	3	Web: 7064, Vid: 274, Tel: 115
VarUnit	0	1.00	FALSE	3	2: 3750, 1: 3689, 3: 14
Stratum	0	1.00	FALSE	50	12: 179, 6: 172, 27: 172, 21: 170
CampaignInterest	1	1.00	FALSE	3	Ver: 3940, Som: 2569, Not: 943
EarlyVote2020	6963	0.07	FALSE	2	Yes: 375, No: 115
VotedPres2016	21	1.00	FALSE	2	Yes: 5810, No: 1622
VotedPres2016_selection	1686	0.77	FALSE	3	Cli: 2911, Tru: 2466, Oth: 390
PartyID	25	1.00	FALSE	7	Str: 1796, Str: 1545, Ind: 881, Ind: 876
TrustGovernment	29	1.00	FALSE	5	Som: 3313, Abo: 2313, Mos: 1016, Nev: 702
TrustPeople	13	1.00	FALSE	5	Mos: 3511, Abo: 2020, Som: 1597, Nev: 264
AgeGroup	294	0.96	FALSE	6	60-: 1436, 70 : 1330, 30-: 1241, 50-: 1200
Education	116	0.98	FALSE	5	Pos: 2514, Bac: 1877, Gra: 1474, Hig: 1160
RaceEth	81	0.99	FALSE	6	Whi: 5420, His: 662, Bla: 650, Asi: 248
Gender	51	0.99	FALSE	2	Fem: 4027, Mal: 3375
Income	517	0.93	FALSE	22	Und: 647, $50: 485, $10: 451, $25: 405
Income7	517	0.93	FALSE	7	$12: 1468, Und: 1076, $20: 1051, $40: 984
VotedPres2020	1053	0.86	FALSE	2	Yes: 6313, No: 87
VotedPres2020_selection	1219	0.84	FALSE	3	Bid: 3509, Tru: 2567, Oth: 158

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
CaseID	0	1.00	336416.23	103653.12	200015.00	225427.00	335416.00	427865.00	535469.00	▇▃▃▆▂
Weight	0	1.00	1.00	1.02	0.01	0.39	0.69	1.21	6.65	▇▂▁▁▁
Age	294	0.96	51.83	17.14	18.00	37.00	53.00	66.00	80.00	▅▇▇▇▇

Residential energy consumption survey

Code

data("recs_2020", package = "srvyrexploR")
recs_2020 %>% 
  select(-matches("^NWEIGHT")) %>% skim()

Data summary
Name	Piped data
Number of rows	18496
Number of columns	39
_______________________
Column type frequency:
character	2
factor	10
logical	2
numeric	25
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
REGIONC	0	1	4	9	0	4	0
STATE_FIPS	0	1	2	2	0	51	0

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
ClimateRegion_BA	0	1.00	FALSE	8	Col: 7116, Mix: 5579, Hot: 2545, Hot: 1577
Urbanicity	0	1.00	FALSE	3	Urb: 12395, Rur: 4081, Urb: 2020
Region	0	1.00	FALSE	4	Sou: 6426, Wes: 4581, Mid: 3832, Nor: 3657
Division	0	1.00	FALSE	10	Sou: 3256, Pac: 2497, Eas: 2014, Mid: 1977
state_postal	0	1.00	FALSE	51	CA: 1152, TX: 1016, NY: 904, FL: 655
state_name	0	1.00	FALSE	51	Cal: 1152, Tex: 1016, New: 904, Flo: 655
HousingUnitType	0	1.00	FALSE	5	Sin: 12319, Apa: 2439, Sin: 1751, Apa: 1013
YearMade	0	1.00	TRUE	9	197: 2817, 200: 2748, Bef: 2721, 199: 2451
HeatingBehavior	751	0.96	FALSE	6	Set: 7806, Man: 4654, Pro: 3310, Tur: 1491
ACBehavior	2325	0.87	FALSE	6	Set: 6738, Man: 3637, Tur: 2746, Pro: 2638

Variable type: logical

skim_variable	n_missing	complete_rate	mean	count
SpaceHeatingUsed	0	1	0.96	TRU: 17745, FAL: 751
ACUsed	0	1	0.87	TRU: 16171, FAL: 2325

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
DOEID	0	1.00	109248.50	5339.48	100001.00	104624.75	109248.50	113872.25	118496.00	▇▇▇▇▇
HDD65	0	1.00	4271.81	2329.46	0.00	2433.75	4396.50	5810.25	17383.00	▆▇▂▁▁
CDD65	0	1.00	1525.52	1143.24	0.00	814.00	1179.00	1805.00	5534.00	▇▆▂▁▁
HDD30YR	0	1.00	4679.41	2338.52	0.00	2897.75	4825.00	6290.00	16071.00	▅▇▃▁▁
CDD30YR	0	1.00	1309.88	988.47	0.00	601.00	1020.00	1703.00	4905.00	▇▅▂▁▁
TOTSQFT_EN	0	1.00	1959.97	1164.97	200.00	1100.00	1700.00	2510.00	15000.00	▇▁▁▁▁
TOTHSQFT	0	1.00	1744.46	1106.33	0.00	1000.00	1520.00	2300.00	15000.00	▇▁▁▁▁
TOTCSQFT	0	1.00	1393.57	1168.31	0.00	460.00	1200.00	2000.00	14600.00	▇▁▁▁▁
WinterTempDay	751	0.96	69.77	3.62	50.00	68.00	70.00	72.00	90.00	▁▁▇▁▁
WinterTempAway	751	0.96	67.45	5.02	50.00	65.00	68.00	70.00	90.00	▁▅▇▁▁
WinterTempNight	751	0.96	68.01	4.66	50.00	65.00	68.00	70.00	90.00	▁▃▇▁▁
SummerTempDay	2325	0.87	72.01	4.74	50.00	70.00	72.00	75.00	90.00	▁▁▇▃▁
SummerTempAway	2325	0.87	73.45	5.68	50.00	70.00	74.00	78.00	90.00	▁▁▇▇▁
SummerTempNight	2325	0.87	71.22	4.80	50.00	68.00	72.00	74.00	90.00	▁▂▇▃▁
BTUEL	0	1.00	37016.17	24265.34	143.32	20205.76	31890.03	48297.98	628155.47	▇▁▁▁▁
DOLLAREL	0	1.00	1424.81	862.08	-889.48	836.49	1257.86	1818.95	15680.18	▇▁▁▁▁
BTUNG	0	1.00	36960.50	46538.32	0.00	0.00	22011.97	62714.20	1134708.69	▇▁▁▁▁
DOLLARNG	0	1.00	396.05	483.41	0.00	0.00	313.86	644.93	8154.98	▇▁▁▁▁
BTULP	0	1.00	3916.95	16741.93	0.00	0.00	0.00	0.00	364215.39	▇▁▁▁▁
DOLLARLP	0	1.00	80.89	333.79	0.00	0.00	0.00	0.00	6621.44	▇▁▁▁▁
BTUFO	0	1.00	5108.60	22698.98	0.00	0.00	0.00	0.00	426268.50	▇▁▁▁▁
DOLLARFO	0	1.00	88.43	395.97	0.00	0.00	0.00	0.00	7003.69	▇▁▁▁▁
BTUWOOD	0	1.00	3596.49	17766.37	0.00	0.00	0.00	0.00	500000.00	▇▁▁▁▁
TOTALBTU	0	1.00	83002.29	53206.00	1182.22	45565.16	74180.11	108535.02	1367548.13	▇▁▁▁▁
TOTALDOL	0	1.00	1990.17	1108.99	-150.51	1258.32	1793.20	2472.00	20043.41	▇▁▁▁▁

Numeric vars

Code

recs_2020 %>% select(where(is.numeric)) %>% colnames()

 [1] "DOEID"           "HDD65"           "CDD65"           "HDD30YR"        
 [5] "CDD30YR"         "TOTSQFT_EN"      "TOTHSQFT"        "TOTCSQFT"       
 [9] "WinterTempDay"   "WinterTempAway"  "WinterTempNight" "SummerTempDay"  
[13] "SummerTempAway"  "SummerTempNight" "NWEIGHT"         "NWEIGHT1"       
[17] "NWEIGHT2"        "NWEIGHT3"        "NWEIGHT4"        "NWEIGHT5"       
[21] "NWEIGHT6"        "NWEIGHT7"        "NWEIGHT8"        "NWEIGHT9"       
[25] "NWEIGHT10"       "NWEIGHT11"       "NWEIGHT12"       "NWEIGHT13"      
[29] "NWEIGHT14"       "NWEIGHT15"       "NWEIGHT16"       "NWEIGHT17"      
[33] "NWEIGHT18"       "NWEIGHT19"       "NWEIGHT20"       "NWEIGHT21"      
[37] "NWEIGHT22"       "NWEIGHT23"       "NWEIGHT24"       "NWEIGHT25"      
[41] "NWEIGHT26"       "NWEIGHT27"       "NWEIGHT28"       "NWEIGHT29"      
[45] "NWEIGHT30"       "NWEIGHT31"       "NWEIGHT32"       "NWEIGHT33"      
[49] "NWEIGHT34"       "NWEIGHT35"       "NWEIGHT36"       "NWEIGHT37"      
[53] "NWEIGHT38"       "NWEIGHT39"       "NWEIGHT40"       "NWEIGHT41"      
[57] "NWEIGHT42"       "NWEIGHT43"       "NWEIGHT44"       "NWEIGHT45"      
[61] "NWEIGHT46"       "NWEIGHT47"       "NWEIGHT48"       "NWEIGHT49"      
[65] "NWEIGHT50"       "NWEIGHT51"       "NWEIGHT52"       "NWEIGHT53"      
[69] "NWEIGHT54"       "NWEIGHT55"       "NWEIGHT56"       "NWEIGHT57"      
[73] "NWEIGHT58"       "NWEIGHT59"       "NWEIGHT60"       "BTUEL"          
[77] "DOLLAREL"        "BTUNG"           "DOLLARNG"        "BTULP"          
[81] "DOLLARLP"        "BTUFO"           "DOLLARFO"        "BTUWOOD"        
[85] "TOTALBTU"        "TOTALDOL"

American national election studies design object

Code

cps_state_in <- getCensus(
  name = "cps/basic/mar",
  vintage = 2020,
  region = "state",
  vars = c(
    "HRMONTH", "HRYEAR4",
    "PRTAGE", "PRCITSHP", "PWSSWGT"
  ),
  key = Sys.getenv("CENSUS_KEY")
)

cps_state <- cps_state_in %>%
  as_tibble() %>%
  mutate(
    across(
      .cols = everything(),
      .fns = as.numeric
    )
  )

Code

cps_narrow_resp <- cps_state %>%
  filter(
    PRTAGE >= 18,
    PRCITSHP %in% c(1:4)
  )

Calculate use population from the narrow data. Weights should add to total pop.

Code

targetpop <- cps_narrow_resp %>% 
  pull(PWSSWGT) %>%
  sum()
scales::comma(targetpop)

[1] "231,034,125"

We can the weight the us election study appropriately

Code

anes_adjwgt <- anes_2020 %>%
  mutate(
    weight = V200010b / sum(V200010b) * targetpop
  )

We then make it conform to a survey design

Code

anes_des <- anes_adjwgt %>%
  as_survey_design(
    weights = weight,
    strata = V200010d,
    ids = V200010c,
    nest = TRUE
  )
anes_des

Stratified 1 - level Cluster Sampling design (with replacement)
With (101) clusters.
Called via srvyr
Sampling variables:
  - ids: V200010c 
  - strata: V200010d 
  - weights: weight 
Data variables: 
  - V200001 (dbl), CaseID (dbl), V200002 (hvn_lbll), InterviewMode (fct),
    V200010b (dbl), Weight (dbl), V200010c (dbl), VarUnit (fct), V200010d
    (dbl), Stratum (fct), V201006 (hvn_lbll), CampaignInterest (fct), V201023
    (hvn_lbll), EarlyVote2020 (fct), V201024 (hvn_lbll), V201025x (hvn_lbll),
    V201028 (hvn_lbll), V201029 (hvn_lbll), V201101 (hvn_lbll), V201102
    (hvn_lbll), VotedPres2016 (fct), V201103 (hvn_lbll),
    VotedPres2016_selection (fct), V201228 (hvn_lbll), V201229 (hvn_lbll),
    V201230 (hvn_lbll), V201231x (hvn_lbll), PartyID (fct), V201233 (hvn_lbll),
    TrustGovernment (fct), V201237 (hvn_lbll), TrustPeople (fct), V201507x
    (hvn_lbll), Age (dbl), AgeGroup (fct), V201510 (hvn_lbll), Education (fct),
    V201546 (hvn_lbll), V201547a (hvn_lbll), V201547b (hvn_lbll), V201547c
    (hvn_lbll), V201547d (hvn_lbll), V201547e (hvn_lbll), V201547z (hvn_lbll),
    V201549x (hvn_lbll), RaceEth (fct), V201600 (hvn_lbll), Gender (fct),
    V201607 (hvn_lbll), V201610 (hvn_lbll), V201611 (hvn_lbll), V201613
    (hvn_lbll), V201615 (hvn_lbll), V201616 (hvn_lbll), V201617x (hvn_lbll),
    Income (fct), Income7 (fct), V202051 (hvn_lbll), V202066 (hvn_lbll),
    V202072 (hvn_lbll), VotedPres2020 (fct), V202073 (hvn_lbll), V202109x
    (hvn_lbll), V202110x (hvn_lbll), VotedPres2020_selection (fct), weight
    (dbl)

Using replicated weights

Code

recs_des <- recs_2020 %>%
  as_survey_rep(
    weights = NWEIGHT,
    repweights = NWEIGHT1:NWEIGHT60,
    type = "JK1",
    scale = 59 / 60,
    mse = TRUE
  )

recs_des

Call: Called via srvyr
Unstratified cluster jacknife (JK1) with 60 replicates and MSE variances.
Sampling variables:
  - repweights: `NWEIGHT1 + NWEIGHT2 + NWEIGHT3 + NWEIGHT4 + NWEIGHT5 +
    NWEIGHT6 + NWEIGHT7 + NWEIGHT8 + NWEIGHT9 + NWEIGHT10 + NWEIGHT11 +
    NWEIGHT12 + NWEIGHT13 + NWEIGHT14 + NWEIGHT15 + NWEIGHT16 + NWEIGHT17 +
    NWEIGHT18 + NWEIGHT19 + NWEIGHT20 + NWEIGHT21 + NWEIGHT22 + NWEIGHT23 +
    NWEIGHT24 + NWEIGHT25 + NWEIGHT26 + NWEIGHT27 + NWEIGHT28 + NWEIGHT29 +
    NWEIGHT30 + NWEIGHT31 + NWEIGHT32 + NWEIGHT33 + NWEIGHT34 + NWEIGHT35 +
    NWEIGHT36 + NWEIGHT37 + NWEIGHT38 + NWEIGHT39 + NWEIGHT40 + NWEIGHT41 +
    NWEIGHT42 + NWEIGHT43 + NWEIGHT44 + NWEIGHT45 + NWEIGHT46 + NWEIGHT47 +
    NWEIGHT48 + NWEIGHT49 + NWEIGHT50 + NWEIGHT51 + NWEIGHT52 + NWEIGHT53 +
    NWEIGHT54 + NWEIGHT55 + NWEIGHT56 + NWEIGHT57 + NWEIGHT58 + NWEIGHT59 +
    NWEIGHT60` 
  - weights: NWEIGHT 
Data variables: 
  - DOEID (dbl), ClimateRegion_BA (fct), Urbanicity (fct), Region (fct),
    REGIONC (chr), Division (fct), STATE_FIPS (chr), state_postal (fct),
    state_name (fct), HDD65 (dbl), CDD65 (dbl), HDD30YR (dbl), CDD30YR (dbl),
    HousingUnitType (fct), YearMade (ord), TOTSQFT_EN (dbl), TOTHSQFT (dbl),
    TOTCSQFT (dbl), SpaceHeatingUsed (lgl), ACUsed (lgl), HeatingBehavior
    (fct), WinterTempDay (dbl), WinterTempAway (dbl), WinterTempNight (dbl),
    ACBehavior (fct), SummerTempDay (dbl), SummerTempAway (dbl),
    SummerTempNight (dbl), NWEIGHT (dbl), NWEIGHT1 (dbl), NWEIGHT2 (dbl),
    NWEIGHT3 (dbl), NWEIGHT4 (dbl), NWEIGHT5 (dbl), NWEIGHT6 (dbl), NWEIGHT7
    (dbl), NWEIGHT8 (dbl), NWEIGHT9 (dbl), NWEIGHT10 (dbl), NWEIGHT11 (dbl),
    NWEIGHT12 (dbl), NWEIGHT13 (dbl), NWEIGHT14 (dbl), NWEIGHT15 (dbl),
    NWEIGHT16 (dbl), NWEIGHT17 (dbl), NWEIGHT18 (dbl), NWEIGHT19 (dbl),
    NWEIGHT20 (dbl), NWEIGHT21 (dbl), NWEIGHT22 (dbl), NWEIGHT23 (dbl),
    NWEIGHT24 (dbl), NWEIGHT25 (dbl), NWEIGHT26 (dbl), NWEIGHT27 (dbl),
    NWEIGHT28 (dbl), NWEIGHT29 (dbl), NWEIGHT30 (dbl), NWEIGHT31 (dbl),
    NWEIGHT32 (dbl), NWEIGHT33 (dbl), NWEIGHT34 (dbl), NWEIGHT35 (dbl),
    NWEIGHT36 (dbl), NWEIGHT37 (dbl), NWEIGHT38 (dbl), NWEIGHT39 (dbl),
    NWEIGHT40 (dbl), NWEIGHT41 (dbl), NWEIGHT42 (dbl), NWEIGHT43 (dbl),
    NWEIGHT44 (dbl), NWEIGHT45 (dbl), NWEIGHT46 (dbl), NWEIGHT47 (dbl),
    NWEIGHT48 (dbl), NWEIGHT49 (dbl), NWEIGHT50 (dbl), NWEIGHT51 (dbl),
    NWEIGHT52 (dbl), NWEIGHT53 (dbl), NWEIGHT54 (dbl), NWEIGHT55 (dbl),
    NWEIGHT56 (dbl), NWEIGHT57 (dbl), NWEIGHT58 (dbl), NWEIGHT59 (dbl),
    NWEIGHT60 (dbl), BTUEL (dbl), DOLLAREL (dbl), BTUNG (dbl), DOLLARNG (dbl),
    BTULP (dbl), DOLLARLP (dbl), BTUFO (dbl), DOLLARFO (dbl), BTUWOOD (dbl),
    TOTALBTU (dbl), TOTALDOL (dbl)

Descriptive statisitcs

Count observations with survey_count() or survey_tally()
Sum variables with survey_total()
Means and proportions: survey_mean(), survey_prop()
Quantiles and medians: survey_quantile(), survey_median()
Correlations: survey_cor()
Ratios: survey_ratio()
Variances and standard deviation: survey_var(), survey_sd()

Counts and cross-tabulations

Code

recs_des %>% survey_tally()

# A tibble: 1 × 2
           n  n_se
       <dbl> <dbl>
1 123529025. 0.148

Estimated counts by subgroups

Code

recs_des %>% survey_count(Region, Division, name = "N")

# A tibble: 10 × 4
   Region    Division                   N         N_se
   <fct>     <fct>                  <dbl>        <dbl>
 1 Northeast New England         5876166  0.0000000137
 2 Northeast Middle Atlantic    16043503  0.0000000487
 3 Midwest   East North Central 18546912  0.000000437 
 4 Midwest   West North Central  8495815  0.0000000177
 5 South     South Atlantic     24843261  0.0000000418
 6 South     East South Central  7380717. 0.114       
 7 South     West South Central 14619094  0.000488    
 8 West      Mountain North      4615844  0.119       
 9 West      Mountain South      4602070  0.0000000492
10 West      Pacific            18505643. 0.00000295

To achieve same result by survey_tally(), you first group

Code

recs_des %>% group_by(Region, Division) %>% survey_tally(name = "N")

# A tibble: 10 × 4
# Groups:   Region [4]
   Region    Division                   N         N_se
   <fct>     <fct>                  <dbl>        <dbl>
 1 Northeast New England         5876166  0.0000000137
 2 Northeast Middle Atlantic    16043503  0.0000000487
 3 Midwest   East North Central 18546912  0.000000437 
 4 Midwest   West North Central  8495815  0.0000000177
 5 South     South Atlantic     24843261  0.0000000418
 6 South     East South Central  7380717. 0.114       
 7 South     West South Central 14619094  0.000488    
 8 West      Mountain North      4615844  0.119       
 9 West      Mountain South      4602070  0.0000000492
10 West      Pacific            18505643. 0.00000295

Totals and sums

To get population count estimate, we leave argument x, empty.

Code

recs_des %>% 
  summarize(
    tot = survey_total()
  )

# A tibble: 1 × 2
         tot tot_se
       <dbl>  <dbl>
1 123529025.  0.148

Overall summation of a continuous variable

Code

recs_des %>%
  summarize(elec_bill = survey_total(DOLLAREL))

# A tibble: 1 × 2
      elec_bill elec_bill_se
          <dbl>        <dbl>
1 170473527909.   664893504.

Summation by groups

Code

recs_des %>% group_by(Region) %>%
  summarize(
  elec_bill = survey_total(DOLLAREL, vartype = "ci")
  
)

# A tibble: 4 × 4
  Region       elec_bill elec_bill_low elec_bill_upp
  <fct>            <dbl>         <dbl>         <dbl>
1 Northeast 29430369947.  28788987554.  30071752341.
2 Midwest   34972544751.  34339576041.  35605513460.
3 South     72496840204.  71534780902.  73458899506.
4 West      33573773008.  32909111702.  34238434313.

You can supply .by argument inside summarize

Code

recs_des  %>%
  summarize(
  elec_bill = survey_total(DOLLAREL, vartype = "ci"),
  .by = Region
)

# A tibble: 4 × 4
  Region       elec_bill elec_bill_low elec_bill_upp
  <fct>            <dbl>         <dbl>         <dbl>
1 Northeast 29430369947.  28788987554.  30071752341.
2 Midwest   34972544751.  34339576041.  35605513460.
3 South     72496840204.  71534780902.  73458899506.
4 West      33573773008.  32909111702.  34238434313.

Means and proportions

Code

recs_des %>% group_by(Region, ACUsed) %>%
  summarise(
  p = survey_prop()
)

When `proportion` is unspecified, `survey_prop()` now defaults to `proportion = TRUE`.
ℹ This should improve confidence interval coverage.
This message is displayed once per session.

# A tibble: 8 × 4
# Groups:   Region [4]
  Region    ACUsed      p    p_se
  <fct>     <lgl>   <dbl>   <dbl>
1 Northeast FALSE  0.110  0.00590
2 Northeast TRUE   0.890  0.00590
3 Midwest   FALSE  0.0666 0.00508
4 Midwest   TRUE   0.933  0.00508
5 South     FALSE  0.0581 0.00278
6 South     TRUE   0.942  0.00278
7 West      FALSE  0.255  0.00759
8 West      TRUE   0.745  0.00759

Joint proportions using interact

Code

recs_des %>% 
  group_by(
    interact(Region, ACUsed)
  ) %>% 
  summarize(
    p = survey_prop() 
  ) %>% mutate(
    p = scales::percent(p)
  )

# A tibble: 8 × 4
  Region    ACUsed p         p_se
  <fct>     <lgl>  <chr>    <dbl>
1 Northeast FALSE  1.96%  0.00105
2 Northeast TRUE   15.79% 0.00105
3 Midwest   FALSE  1.46%  0.00111
4 Midwest   TRUE   20.43% 0.00111
5 South     FALSE  2.20%  0.00106
6 South     TRUE   35.72% 0.00106
7 West      FALSE  5.73%  0.00170
8 West      TRUE   16.72% 0.00170

Overall mean

Code

recs_des %>% summarize(
  elec_bill = survey_mean(DOLLAREL, vartype = c("se", "ci"))
)

# A tibble: 1 × 4
  elec_bill elec_bill_se elec_bill_low elec_bill_upp
      <dbl>        <dbl>         <dbl>         <dbl>
1     1380.         5.38         1369.         1391.

Quantiles and Medians

Overall quantiles

Code

recs_des %>%
  summarize(
    elec_bill = survey_quantile(DOLLAREL, quantiles = c(.25, .5, .75))
  )

# A tibble: 1 × 6
  elec_bill_q25 elec_bill_q50 elec_bill_q75 elec_bill_q25_se elec_bill_q50_se
          <dbl>         <dbl>         <dbl>            <dbl>            <dbl>
1          795.         1215.         1770.             5.69             6.33
# ℹ 1 more variable: elec_bill_q75_se <dbl>

Overall median

Code

recs_des %>%
  summarize(
    elec_bill = survey_median(DOLLAREL)
  )

# A tibble: 1 × 2
  elec_bill elec_bill_se
      <dbl>        <dbl>
1     1215.         6.33

Correlations

Code

recs_des %>%
  summarize(SQFT_Elec_Corr = survey_corr(TOTSQFT_EN, BTUEL))

Warning: There was 1 warning in `dplyr::summarise()`.
ℹ In argument: `SQFT_Elec_Corr = survey_corr(TOTSQFT_EN, BTUEL)`.
Caused by warning in `sweep()`:
! length(STATS) or dim(STATS) do not match dim(x)[MARGIN]

# A tibble: 1 × 2
  SQFT_Elec_Corr SQFT_Elec_Corr_se
           <dbl>             <dbl>
1          0.417           0.00689

Design effects

Design effect measures the precision of an estimate under a particular sampling design relative to a simple random sampling design. If <1, the design is statistically more efficient than SRS. It is used in the calculation of effective sample size - the sample size needed if we were to employ SRS

\[n_{eff} = \frac{n}{D_{eff}}\]

Code

recs_des %>%
  summarize(
    across(
          c(BTUEL, BTUNG, BTULP, BTUFO, BTUWOOD),
          ~ survey_mean(.x, deff = TRUE, vartype = NULL)
    )
  ) %>% select(ends_with("deff"))

# A tibble: 1 × 5
  BTUEL_deff BTUNG_deff BTULP_deff BTUFO_deff BTUWOOD_deff
       <dbl>      <dbl>      <dbl>      <dbl>        <dbl>
1      0.597      0.938       1.21      0.720         1.10

Creating cascading rows - `cascade`

Code

recs_des %>%
  cascade(DOLLAREL_mn = survey_mean(DOLLAREL))

# A tibble: 1 × 2
  DOLLAREL_mn DOLLAREL_mn_se
        <dbl>          <dbl>
1       1380.           5.38

Group

Code

recs_des %>%
  group_by(Region) %>%
  cascade(DOLLAREL_mn = survey_mean(DOLLAREL), .fill = "National")

# A tibble: 5 × 3
  Region    DOLLAREL_mn DOLLAREL_mn_se
  <fct>           <dbl>          <dbl>
1 Northeast       1343.          14.6 
2 Midwest         1293.          11.7 
3 South           1548.          10.3 
4 West            1211.          12.0 
5 National        1380.           5.38

Or through a long dplyr pipeline

Code

recs_des %>% 
  summarize(
    cons_mean = survey_mean(DOLLAREL),
    .by = Region
  ) %>% bind_rows(
    ., recs_des %>%
      summarize(
        cons_mean = survey_mean(DOLLAREL)
      ) %>% mutate(
        ., Region = "National"
      ) %>% select(Region, everything())
  )

# A tibble: 5 × 3
  Region    cons_mean cons_mean_se
  <chr>         <dbl>        <dbl>
1 Northeast     1343.        14.6 
2 Midwest       1293.        11.7 
3 South         1548.        10.3 
4 West          1211.        12.0 
5 National      1380.         5.38

Summarize across multiple variables

Code

recs_des %>% 
  summarize(
    across(
      starts_with("BTU"),
      list(
        Total = ~ survey_total(.x, vartype = "cv"),
        Mean = ~ survey_mean(.x, vartype = "cv")
      ),
      .unpack = "{outer}.{inner}"
    )
  )

# A tibble: 1 × 20
  BTUEL_Total.coef BTUEL_Total._cv BTUEL_Mean.coef BTUEL_Mean._cv
             <dbl>           <dbl>           <dbl>          <dbl>
1    4453284510065         0.00377          36051.        0.00377
# ℹ 16 more variables: BTUNG_Total.coef <dbl>, BTUNG_Total._cv <dbl>,
#   BTUNG_Mean.coef <dbl>, BTUNG_Mean._cv <dbl>, BTULP_Total.coef <dbl>,
#   BTULP_Total._cv <dbl>, BTULP_Mean.coef <dbl>, BTULP_Mean._cv <dbl>,
#   BTUFO_Total.coef <dbl>, BTUFO_Total._cv <dbl>, BTUFO_Mean.coef <dbl>,
#   BTUFO_Mean._cv <dbl>, BTUWOOD_Total.coef <dbl>, BTUWOOD_Total._cv <dbl>,
#   BTUWOOD_Mean.coef <dbl>, BTUWOOD_Mean._cv <dbl>

Map

Code

calc <- function(df, var) {
  df %>%
    drop_na(!!sym(var)) %>%
    group_by(!!sym(var)) %>%
    summarize(p = survey_prop()) %>%
    mutate(Variable = var) %>%
    rename(Answer := !!sym(var)) %>%
    mutate(p = scales::percent(p)) %>%
    select(Variable, everything())
}

v <- c("TrustGovernment", "TrustPeople")
map2(
  .x = list(anes_des),
  .y = v,
  .f = ~calc(df = .x, var = .y)
) %>% bind_rows()

# A tibble: 10 × 4
   Variable        Answer              p        p_se
   <chr>           <fct>               <chr>   <dbl>
 1 TrustGovernment Always              1.6%  0.00204
 2 TrustGovernment Most of the time    13.2% 0.00553
 3 TrustGovernment About half the time 30.9% 0.00829
 4 TrustGovernment Some of the time    43.4% 0.00855
 5 TrustGovernment Never               11.0% 0.00566
 6 TrustPeople     Always              0.8%  0.00164
 7 TrustPeople     Most of the time    41.4% 0.00857
 8 TrustPeople     About half the time 28.2% 0.00776
 9 TrustPeople     Some of the time    24.5% 0.00670
10 TrustPeople     Never               5.0%  0.00422

Statistical testing

Comparison of means and proportions: svyttest()
goodness-of-fit test: svygofchisq()
Test of independence: svychisq()
Test of homogeneity: svychisq()

Comparison of Means and Proportions

one-sample t-test: compare to 0 - var ~ 0 or to a different value var - val = 0
two-sample t-test:
- unpaired: two level grouping variable - var ~ groupvar or three level grouping variable - var ~ groupvar == level
- paired: var_1 - var_2 = 0

Example: one-sample t-test for mean

Code

recs_des %>%
  svyttest(
    formula = SummerTempNight - 68 ~ 0,
    design =.,
    na.rm = TRUE
  )


    Design-based one-sample t-test

data:  SummerTempNight - 68 ~ 0
t = 84.788, df = 58, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 3.287816 3.446810
sample estimates:
    mean 
3.367313

get mean summer temps at night

Code

recs_des %>% 
  summarize(m = survey_mean(SummerTempNight, na.rm = TRUE))

# A tibble: 1 × 2
      m   m_se
  <dbl>  <dbl>
1  71.4 0.0397

Example: one-sample t-test for proportion

Code

recs_des %>% summarize(
  p = survey_prop(),
  .by = ACUsed
)

# A tibble: 2 × 3
  ACUsed     p    p_se
  <lgl>  <dbl>   <dbl>
1 FALSE  0.113 0.00306
2 TRUE   0.887 0.00306

Code

recs_des %>%
  svyttest(
    formula = (ACUsed == TRUE) - 0.90 ~ 0,
    design = .,
    na.rm = TRUE
  ) %>% tidy() %>%
  mutate(p.value = pretty_p_value(p.value)) %>%
  gt() %>%
  fmt_number()

Table 1: One-sample t-test for estimate of US housholds use of A/C in their homes differing from 90%

estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative
−0.01	−4.40	<0.0001	58.00	−0.02	−0.01	Design-based one-sample t-test	two.sided

Example: paired two-sample t-test

Code

recs_des %>%
  svyttest(
    formula = SummerTempNight - WinterTempNight ~ 0,
    design = .,
    na.rm = TRUE
  ) %>% tidy() %>%
  mutate(
    p.value = pretty_p_value(p.value)
  ) %>%
  gt() %>%
  fmt_number()

Table 2: Paired two-sample t-test

estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative
2.85	50.83	<0.0001	58.00	2.74	2.96	Design-based one-sample t-test	two.sided

Chi-squared tests

See https://tidy-survey-r.github.io/tidy-survey-book/c06-statistical-testing.html

Example: goodness of fit test

Code

anes_des_educ <- anes_des %>%
  mutate(
    Education2 = fct_collapse(Education,
      "Bachelor or Higher" = c("Bachelor's", "Graduate")
    )
  )

anes_des_educ %>%
  drop_na(Education2) %>%
  group_by(Education2) %>% 
  summarize(p = survey_mean())

# A tibble: 4 × 3
  Education2              p    p_se
  <fct>               <dbl>   <dbl>
1 Less than HS       0.0805 0.00568
2 High school        0.277  0.0102 
3 Post HS            0.290  0.00713
4 Bachelor or Higher 0.352  0.00732

Code

anes_des_educ %>%
  svygofchisq(
    formula = ~Education2,
    design = .,
    p = c(0.11, 0.27, 0.28, 0.35),
    na.rm = TRUE
  ) %>% tidy()

Multiple parameters; naming those columns scale and df.

# A tibble: 1 × 5
    scale    df statistic  p.value method                                       
    <dbl> <dbl>     <dbl>    <dbl> <chr>                                        
1 113213.  2.26  1922280. 0.000292 Design-based chi-squared test for given prob…

Code

ex1_table <- anes_des_educ %>%
  drop_na(Education2) %>%
  group_by(Education2) %>%
  summarize(Observed = survey_mean(vartype = "ci")) %>%
  rename(Education = Education2) %>%
  mutate(Expected = c(0.11, 0.27, 0.29, 0.33)) %>%
  select(Education, Expected, everything())

ex1_table

# A tibble: 4 × 5
  Education          Expected Observed Observed_low Observed_upp
  <fct>                 <dbl>    <dbl>        <dbl>        <dbl>
1 Less than HS           0.11   0.0805       0.0691       0.0919
2 High school            0.27   0.277        0.257        0.298 
3 Post HS                0.29   0.290        0.276        0.305 
4 Bachelor or Higher     0.33   0.352        0.337        0.367

Code

ex1_table %>%
  pivot_longer(
    cols = c("Expected", "Observed"),
    names_to = "Names",
    values_to = "Proportion"
  ) %>%
  mutate(
    Observed_low = if_else(Names == "Observed", Observed_low, NA_real_),
    Observed_upp = if_else(Names == "Observed", Observed_upp, NA_real_),
    Names = if_else(Names == "Observed",
      "ANES (observed)", "ACS (expected)"
    )
  ) %>%
  ggplot(aes(x = Education, y = Proportion, color = Names)) +
  geom_point(alpha = 0.75, size = 2) +
  geom_errorbar(aes(ymin = Observed_low, ymax = Observed_upp),
    width = 0.25
  ) +
  theme_bw() +
  # scale_color_manual(name = "Type", values = book_colors[c(4, 1)]) +
  theme(legend.position = "bottom", legend.title = element_blank())

Sampling Designs and replicate weights

SRS

SRS without replacement

When we sample without replacement, we use finite population correction in the estimation of standard error of estimates

\[fpc = \sqrt{\frac{N-n}{N-1}} \simeq \sqrt{1 - \frac{n}{N}}\]

Code

set.seed(as.numeric(Sys.time()))

data("api", package = "survey")

apisrs_slim <-
  apisrs %>%
  as_tibble() %>%
  arrange(dnum, snum) %>%
  select(cds, dnum, snum, dname, sname, fpc, pw)

apisrs_slim

# A tibble: 200 × 7
   cds             dnum  snum dname                   sname            fpc    pw
   <chr>          <int> <dbl> <chr>                   <chr>          <dbl> <dbl>
 1 19642126061220     1  1121 ABC Unified             Haskell (Plin…  6194  31.0
 2 19642126066716     1  1124 ABC Unified             Stowers (Ceci…  6194  31.0
 3 36675876035174     5  3895 Adelanto Elementary     Adelanto Elem…  6194  31.0
 4 33669776031512    19  3347 Alvord Unified          Arlanza Eleme…  6194  31.0
 5 33669776031595    19  3352 Alvord Unified          Wells Interme…  6194  31.0
 6 31667876031033    39  3271 Auburn Union Elementary Cain (E.V.) M…  6194  31.0
 7 19642876011407    42  1169 Baldwin Park Unified    Deanza Elemen…  6194  31.0
 8 19642876011464    42  1175 Baldwin Park Unified    Heath (Margar…  6194  31.0
 9 19642956011589    48  1187 Bassett Unified         Erwin (Thomas…  6194  31.0
10 41688586043392    49  4948 Bayshore Elementary     Bayshore Elem…  6194  31.0
# ℹ 190 more rows

Create a tbl_survey object

Code

apisrs_des <- apisrs_slim %>%
  as_survey_design(
    weights = pw,
    fpc = fpc
  )
apisrs_des

Independent Sampling design
Called via srvyr
Sampling variables:
  - ids: `1` 
  - fpc: fpc 
  - weights: pw 
Data variables: 
  - cds (chr), dnum (int), snum (dbl), dname (chr), sname (chr), fpc (dbl), pw
    (dbl)

SRS with replacement

Code

set.seed(as.numeric(Sys.time()))

apisrswr <- apipop %>%
  as_tibble() %>%
  slice_sample(n = 200, replace = TRUE) %>%
  select(cds, dnum, snum, dname, sname) %>%
  mutate(weight = nrow(apipop) / 200)

head(apisrswr)

# A tibble: 6 × 6
  cds             dnum  snum dname                     sname              weight
  <chr>          <int> <dbl> <chr>                     <chr>               <dbl>
1 01611926001051   294   100 Hayward Unified           Lorin A. Eden Ele…   31.0
2 10621666006266   253   645 Fresno Unified            Gibson Elementary    31.0
3 36678506036669   590  4101 Rialto Unified            Morgan Elementary    31.0
4 35674706071088   304  3886 Hollister School District Rancho San Justo …   31.0
5 16639821635606   378  1064 Lemoore Union High        Lemoore High         31.0
6 37683616040372   668  4610 Santee Elementary         Pepper Drive Elem…   31.0

Statified Sampling

Population can be grouped into homogeneous units (strata). Random samples are then drawn from each of the units.

If $\hat{y}_h$ is the sample mean for stratum $h$, $N_h$ the population size of stratum $h$, and $H$ the total number of strata, then estimates of the mean

\[\hat{y} = \frac{1}{N}\sum_{h=1}^H N_h \hat{y}_h\] and

\[se(\hat{y}) = \sqrt{\frac{1}{N^2}\sum_{h=1}^H N^2_h \frac{s^2_h}{n_h}(1 - \frac{n_h}{N_h})}\] where

\[s^2_h = \frac{1}{n_h - 1} \sum_{i=1}^{n_h}(y_{i,h} - \hat{y_h})^2\] and proportion

\[\hat{p} = \frac{1}{N} \sum_{h=1}^H N_h \hat{p}_h\]

\[se(\hat{p}) = \frac{1}{N}\sqrt{\sum_{h=1}^H N^2_h\frac{\hat{p_h}(1- \hat{p_h})}{n_h - 1}(1 - \frac{n_h}{N_h})}\]

Clustered sampling

Applies if a population can be divided into mutually exclusive subgroups (clusters or primary sampling units (PSU)).A random selection of PSUs is sampled followed by another level of sampling within the chosen clusters. Suppose $a$ clusters are sampled from a population of A clusters via SRS. Within each sampled cluster $i$ there are $B_i$ units in the population of which $b_i$ units are sampled using SRS. if $\hat{y}_i$ is the $i$th cluster mean, then the population mean is estimated by

\[\hat{y} = \frac{\sum_{i=1}^aB_i\hat{y_i}}{\sum_{i=1}^a B_i}\]

\[se(\hat{y}) = \frac{1}{\hat{N}}\sqrt{(1 - \frac{a}{A})\frac{s^2_a}{a} + \frac{A}{a}\sum_{i=1}^a(1 - \frac{b_i}{B_i}) \frac{s_i^2}{b_i}}\] where $\hat{N}$ is the estimated population size, $s^2_a$ is the between-cluster variance, and $s^2_i$ is the within-cluster variance.

Between-cluster variance $s^2_a$ is:

\[s_a^2 = \frac{1}{a - 1} \sum_{i = 1}^a (\hat{y_i} - \frac{\sum_{i=1}^a \hat{y_i}}{a})^2\] where $\frac{\sum_{i=1}^a \hat{y_i}}{a} = \bar{y}$ is the grand mean.

The within-cluster variance ($s^2_i$) is estimated through

\[s^2_i = \frac{1}{a(b_i - 1)} \sum_{j = 1}^{b_i}(y_{ij} - \hat{y_i})\] where $y_{ij}$ is the outcome for sampled unit $j$ from cluster $i$

Citation

BibTeX citation:

@online{okola2026,
  author = {Okola, Basil},
  title = {Surveys: {Design} and {Analysis}},
  date = {2026-01-15},
  url = {https://bokola.github.io/posts/2026-01-15-survey-design/},
  langid = {en}
}

For attribution, please cite this work as:

Okola, Basil. 2026. “Surveys: Design and Analysis.” January 15, 2026. https://bokola.github.io/posts/2026-01-15-survey-design/.

Analysis

Data

Descriptive statisitcs

Counts and cross-tabulations

Totals and sums

Means and proportions

Quantiles and Medians

Design effects

Creating cascading rows - cascade

Summarize across multiple variables

Map

Statistical testing

Comparison of Means and Proportions

Chi-squared tests

Sampling Designs and replicate weights

SRS

SRS without replacement

SRS with replacement

Statified Sampling

Clustered sampling

Citation

Creating cascading rows - `cascade`