## Survey Data Analysis with R

### Why do we need survey data analysis software?

Regular procedures in statistical software (that is not designed for survey data) analyzes data as if the data were collected using simple random sampling. For experimental and quasi-experimental designs, this is exactly what we want. However, very few surveys use a simple random sample to collect data. Not only is it nearly impossible to do so, but it is not as efficient (either financially and statistically) as other sampling methods. When any sampling method other than simple random sampling is used, we usually need to use survey data analysis software to take into account the differences between the design that was used to collect the data and simple random sampling. This is because the sampling design affects both the calculation of the point estimates and the standard errors of those estimates. If you ignore the sampling design, e.g., if you assume simple random sampling when another type of sampling design was used, both the point estimates and their standard errors will likely be calculated incorrectly. The sampling weight will affect the calculation of the point estimate, and the stratification and/or clustering will affect the calculation of the standard errors. Ignoring the clustering will likely lead to standard errors that are underestimated, possibly leading to results that seem to be statistically significant, when in fact, they are not. The difference in point estimates and standard errors obtained using non-survey software and survey software with the design properly specified will vary from data set to data set, and even between analyses using the same data set. While it may be possible to get reasonably accurate results using non-survey software, there is no practical way to know beforehand how far off the results from non-survey software will be.

### Sampling designs

Most people do not conduct their own surveys. Rather, they use survey data that some agency or company collected and made available to the public. The documentation must be read carefully to find out what kind of sampling design was used to collect the data. This is very important because many of the estimates and standard errors are calculated differently for the different sampling designs. Hence, if you mis-specify the sampling design, the point estimates and standard errors will likely be wrong.

Below are some common features of many sampling designs.

__Sampling weights__: There are several types of weights that can be associated with a survey. Perhaps the most common is the sampling weight. A sampling weight is a probability weight that has had one or more adjustments made to it. Both a sampling weight and a probability weight are used to weight the sample back to the population from which the sample was drawn. By definition, a probability weight is the inverse of the probability of being included in the sample due to the sampling design (except for a certainty PSU, see below). The probability weight is calculated as N/n, where N = the number of elements in the population and n = the number of elements in the sample. For example, if a population has 10 elements and 3 are sampled at random with replacement, then the probability weight would be 10/3 = 3.33. In a two-stage design, the probability weight is calculated as f_{1}f_{2}, which means that the inverse of the sampling fraction for the first stage is multiplied by the inverse of the sampling fraction for the second stage. Under many sampling plans, the sum of the probability weights will equal the population total.

While many textbooks will end their discussion of probability weights here, this definition does not fully describe the sampling weights that are included with actual survey data sets. Rather, the sampling weight, which is sometimes called a “final weight,” starts with the inverse of the sampling fraction, but then incorporates several other values, such as corrections for unit non-response, errors in the sampling frame (sometimes called non-coverage), calibration and trimming. Because these other values are included in the probability weight that is included with the data set, it is often inadvisable to modify the sampling weights, such as trying to standardize them for a particular variable, e.g., age.

__Strata__: Stratification is a method of breaking up the population into different groups, often by demographic variables such as gender, race or SES. Each element in the population must belong to one, and only one, strata. Once the strata have been defined, samples are taken from each stratum as if it were independent of all of the other strata. For example, if a sample is to be stratified on gender, men and women would be sampled independently of one another. This means that the probability weights for men will likely be different from the probability weights for the women. In most cases, you need to have two or more PSUs in each stratum. The purpose of stratification is to reduce the standard error of the estimates, and stratification works most effectively when the variance of the dependent variable is smaller within the strata than in the sample as a whole.

__PSU__: This is the **p**rimary **s**ampling **u**nit. This is the first unit that is sampled in the design. For example, school districts from California may be sampled and then schools within districts may be sampled. The school district would be the PSU. If states from the US were sampled, and then school districts from within each state, and then schools from within each district, then states would be the PSU. One does not need to use the same sampling method at all levels of sampling. For example, probability-proportional-to-size sampling may be used at level 1 (to select states), while cluster sampling is used at level 2 (to select school districts). In the case of a simple random sample, the PSUs and the elementary units are the same. In general, accounting for the clustering in the data (i.e., using the PSUs), will increase the standard errors of the point estimates. Conversely, ignoring the PSUs will tend to yield standard errors that are too small, leading to false positives when doing significance tests.

__FPC__: This is the **f**inite **p**opulation **c**orrection. This is used when the sampling fraction (the number of elements or respondents sampled relative to the population) becomes large. The FPC is used in the calculation of the standard error of the estimate. If the value of the FPC is close to 1, it will have little impact and can be safely ignored. In some survey data analysis programs, such as SUDAAN, this information will be needed if you specify that the data were collected without replacement (see below for a definition of “without replacement”). The formula for calculating the FPC is ((N-n)/(N-1))^{1/2}, where N is the number of elements in the population and n is the number of elements in the sample. To see the impact of the FPC for samples of various proportions, suppose that you had a population of 10,000 elements.

Sample size (n) FPC

1 1.0000

10 .9995

100 .9950

500 .9747

1000 .9487

5000 .7071

9000 .3162

__Replicate weights__: Replicate weights are a series of weight variables that are used to correct the standard errors for the sampling plan. They serve the same function as the PSU and strata variables (which are used a Taylor series linearization) to correct the standard errors of the estimates for the sampling design. Many public use data sets are now being released with replicate weights instead of PSUs and strata in an effort to more securely protect the identity of the respondents. In theory, the same standard errors will be obtained using either the PSU and strata or the replicate weights. There are different ways of creating replicate weights; the method used is determined by the sampling plan. The most common are balanced repeated and jackknife replicate weights. You will need to read the documentation for the survey data set carefully to learn what type of replicate weight is included in the data set; specifying the wrong type of replicate weight will likely lead to incorrect standard errors. For more information on replicate weights, please see Stata Library: Replicate Weights. Several statistical packages, including Stata, SAS, R, Mplus, SUDAAN and WesVar, allow the use of replicate weights. Another good source of information on replicate weights is *Applied Survey Data Analysis, Second Edition* by Steven G. Heeringa, Brady T. West and Patricia A. Berglund (2017, CRC Press).

### Consequences of not using the design elements

Sampling design elements include the sampling weights, post-stratification weights (if provided), PSUs, strata, and replicate weights. Rarely are all of these elements included in a single public-use data set. However, ignoring the design elements that are included can often lead to inaccurate point estimates and/or inaccurate standard errors.

### Sampling with and without replacement

Most samples collected in the real world are collected “without replacement”. This means that once a respondent has been selected to be in the sample and has participated in the survey, that particular respondent cannot be selected again to be in the sample. Many of the calculations change depending on if a sample is collected with or without replacement. Hence, programs like SUDAAN request that you specify if a survey sampling design was implemented with our without replacement, and an FPC is used if sampling without replacement is used, even if the value of the FPC is very close to one.

### Examples

For the examples in this workshop, we will use the data set from NHANES 2011-2012. The data set and documentation can be downloaded from the NHANES web site. The data files can be downloaded as SAS.xpt files. The R script used for the workshop is here.

### Reading the documentation

The first step in analyzing any survey data set is to read the documentation. With many of the public use data sets, the documentation can be quite extensive and sometimes even intimidating. Instead of trying to read the documentation “cover to cover”, there are some parts you will want to focus on. First, read the Introduction. This is usually an “easy read” and will orient you to the survey. There is usually a section or chapter called something like “Sample Design and Analysis Guidelines”, “Variance Estimation”, etc. This is the part that tells you about the design elements included with the survey and how to use them. Some even give example code. If multiple sampling weights have been included in the data set, there will be some instruction about when to use which one. If there is a section or chapter on missing data or imputation, please read that. This will tell you how missing data were handled. You should also read any documentation regarding the specific variables that you intend to use. As we will see little later on, we will need to look at the documentation to get the value labels for the variables. This is especially important because some of the values are actually missing data codes, and you need to do something so that R doesn’t treat those as valid values (or you will get some very “interesting” means, totals, etc.).

### Downloading and installing the packages

You should use the following **install.packages** commands only once. The commands should be left in your R script file but commented out (by placing a “#” symbol before each).

install.packages("haven")install.packages("survey")install.packages("jtools")install.packages("remotes")remotes::install_github("carlganz/svrepmisc")

After the packages are downloaded, they need to be loaded. This needs to be done at the beginning of each R session.

library("haven")library("survey")library("jtools")library("remotes")library("svrepmisc")

### Reading the data into R

The data are distributed from the NHANES website as SAS.xpt files. You can use the “foreign” package to read the data into R. The command will look like this:

nhanes2012 <- read_dta("D:/data/Seminars/Applied Survey Data Analysis R/nhanes2012_reduced.dta")

### The variables

We will use about a dozen different variables in the examples in this workshop. Below is a brief summary of them. Some of the variables have been recoded to be binary variables (values of 2 recoded to a value of 0). The count of missing observations includes values truly missing as well as refused and don’t know.

**ridageyr** – Age in years at exam – recoded; range of values: 0 – 79 are actual values, 80 = 80+ years of age

**pad630** – How much time do you spend doing moderate-intensity activities on a type work day?; range of values: 10-960, 7053 missing observations

**hsq496** – During the past 30 days, for about how many days have you felt worried, tense or anxious?; range of values: 0-30; 3073 missing observations

**female** – Recode of the variable riagendr; 0 = male, 1 = female; no missing observations

**dmdborn4** – Country of birth; 1 = born in the United States, 0 = otherwise; 5 missing observations

**dmdmartl** – Marital status; 1 = married, 2 = widowed, 3 = divorced, 4 = separated, 5 = never married, 6 = living with partner; 4203 missing observations

**dmdeduc2** – Education level of adults aged 20+ years; 1 = less than 9th grade, 2 = 9-11th grade, 3 = high school graduate, GED or equivalent, 4 = some college or AA degree, 5 = college graduate or above; 4201 missing observations

**pad675** – How much time do you spend doing moderate-intensity sports, fitness, or recreation activities on a typical day?; range of values: 10-600; 6220 missing observations

**hsq571** – During the past 12 months, have you donated blood?; 0 = no, 1 = yes; 3673 missing observations

**pad680** – How much time do you usually spend sitting on a typical day?; range of values: 0-1380; 2365 missing observations

**paq665** – Do you do any moderate-intensity sports, fitness or recreational activities that cause a small increase in breathing or heart rate at least 10 minutes continually?; 0 = no, 1 = yes; 2329 missing observations

**hsd010** – Would you say that your general health is…; 1 = excellent, 2 = very good, 3 = good, 4 = fair, 5 = poor; 3064 missing observations

### The svydesign function

Before we can start our analyses, we need to use the **svydesign** function from the “survey” package written by Thomas Lumley. The **svydesign** function tells R about the design elements in the survey. Once this command has been issued, all you need to do for your analyses is use the object that contains this information in each command. Because the 2011-2012 NHANES data were released with a sampling weight (**wtint2yr**), a PSU variable (**sdmvpsu**) and a strata variable (**sdmvstra**), we will use these our **svydesign** function. The **svydesign** function looks like this:

# nhanes2012 <- read_dta("D:/data/Seminars/Applied Survey Data Analysis R/nhanes2012_reduced.dta") nhc <- svydesign(id=~sdmvpsu, weights=~wtint2yr,strata=~sdmvstra, nest=TRUE, survey.lonely.psu = "adjust", data=nhanes2012) nhc

Stratified 1 - level Cluster Sampling design (with replacement) With (31) clusters.

svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012)

We can get additional information about the sample, such as the number of PSUs per strata, by using the **summary** function.

summary(nhc)Stratified 1 - level Cluster Sampling design (with replacement) With (31) clusters. svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Probabilities: Min. 1st Qu. Median Mean 3rd Qu. Max. 4.541e-06 2.866e-05 5.526e-05 6.372e-05 8.809e-05 3.011e-04 Stratum Sizes: 90 91 92 93 94 95 96 97 98 99 100 101 102 103 obs 862 998 875 602 688 722 676 608 708 682 700 715 624 296 design.PSU 3 3 3 2 2 2 2 2 2 2 2 2 2 2 actual.PSU 3 3 3 2 2 2 2 2 2 2 2 2 2 2 Data variables: [1] "dmdborn4" "dmdeduc2" "dmdmartl" "dmdyrsus" "female" "hsd010" "hsq496" [8] "hsq571" "pad615" "pad630" "pad660" "pad675" "pad680" "paq610" [15] "paq625" "paq640" "paq655" "paq665" "paq670" "paq710" "paq715" [22] "ridageyr" "sdmvpsu" "sdmvstra" "seqn" "wtint2yr"

### Descriptive statistics with continuous variables

We will start with something simple: calculating the mean of a continuous variable. In this example, we use the variable **ridageyr**, which is the age of the respondent. Please note that the documentation for the **svymean** function, as well as other functions that provide descriptive statistics, is found in the section of the documentation called **surveysummary**.

svymean(~ridageyr, nhc)mean SE ridageyr 37.185 0.6965

We can also get the standard deviation of the **age** variable. We use the function **svysd**, which is found in the **jtools** package.

svysd(~ridageyr,design = nhc, na = TRUE)std. dev. ridageyr 22.37

When there are missing data for a variable, the **na = TRUE** argument is needed.

svymean(~pad630, nhc, na = TRUE)mean SE pad630 139.89 5.5791

Here is another example.

svymean(~hsq496, nhc, na = TRUE)mean SE hsq496 5.3839 0.19

The means of more than one variable can be obtained by placing “+” between the variables. Notice that a listwise deletion has been done, so that the means in output below are different from the means shown above.

svymean(~ridageyr+pad630+hsq496, nhc, na = TRUE)mean SE ridageyr 41.9779 0.8212 pad630 139.7228 5.7137 hsq496 5.3739 0.2781

In the examples below, the mean, standard deviation, and variance are obtained.

svymean(~pad680, nhc, na = TRUE)mean SE pad680 391.3 5.986

svysd(~pad680,design = nhc, na = TRUE)std. dev. pad680 200.846

svyvar(~pad680, design = nhc, na = TRUE)variance SE pad680 40339 827.77

The **cv** function is used to get the coefficient of variation. The coefficient of variation is the ratio of the standard error to the mean, multiplied by 100%. It is an indication of the variability relative to the mean in the population and is not affected by the unit of measurement of the variable.

cv(svymean(~pad680,design = nhc, na = TRUE))pad680 pad680 0.01529783

Below is an example of how to obtain the design effect (Deff). The Deff is a ratio of two variances. In the numerator we have the variance estimate from the current sample (including all of its design elements), and in the denominator we have the variance from a hypothetical sample of the same size drawn as an SRS. In other words, the Deff tells you how efficient your sample is compared to an SRS of equal size. If the Deff is less than 1, your sample is more efficient than SRS; usually the Deff is greater than 1. In the example below, the Deff is 5.9887. This means that a sample drawn using the current sampling plan needs to be six times the size needed if the sample was collected via an SRS. The Deff is specific to a variable, so some variables may have larger or smaller Deffs.

svymean(~pad680, nhc, na = TRUE, deff = TRUE)mean SE DEff pad680 391.300 5.986 5.9887

Quantiles are a useful descriptive statistic for continuous variables, particularly variables that are not normally distributed.

svyquantile(~hsq496, design = nhc, na = TRUE, c(.25,.5,.75),ci=TRUE)$quantiles 0.25 0.5 0.75 hsq496 0 1 5 $CIs , , hsq496 0.25 0.5 0.75 (lower 0 1 5 upper) 0 2 7

Let’s get some descriptive statistics for binary variables. We have some choices here. Let’s start by getting the mean of the variable **female**, which is coded 1 for females and 0 for males. Taking the mean of a variable that is coded 0/1 gives the proportion of 1s, so the mean of this variable is the estimated proportion of the population that is female.

# means and proportions for binary variablessvymean(~female, nhc)mean SE female 0.51195 0.0064

We can use the **confint** function to get the confidence interval around this mean.

confint(svymean(~female, nhc))2.5 % 97.5 % female 0.4993307 0.5245742

However, comment on page 70of the documentation for the **survey** package, we should use **svyciprop** rather than **confint**. There are several options that can be supplied for the **method** argument. Please see pages 70-71 of the documentation. The **likelihood** option uses the (Rao-Scott) scaled chi-squared distribution for the log likelihood from a binomial distribution.

svyciprop(~I(female==1), nhc, method="likelihood")2.5% 97.5% I(female == 1) 0.512 0.498 0.53

# li is short for likelihoodsvyciprop(~I(female==0), nhc, method="li")2.5% 97.5% I(female == 0) 0.488 0.474 0.5

The **logit** option fits a logistic regression model and computes a Wald-type interval on the log-odds scale, which is then transformed to the probability scale.

svyciprop(~I(female==1), nhc, method="logit")2.5% 97.5% I(female == 1) 0.512 0.498 0.53

The **xlogit** option uses a logit transformation of the mean and then back-transforms to the probability scale. This appears to be the method used by SUDAAN and SPSS COMPLEX SAMPLES.

svyciprop(~I(female==1), nhc, method="xlogit")2.5% 97.5% I(female == 1) 0.512 0.498 0.53

As stated on page 71 of the documentation, the use of the mean option (shortened in the code below to me), reproduces the results given by Stata’s **svy: mean** command.

svyciprop(~I(female==1), nhc, method="me", df=degf(nhc))2.5% 97.5% I(female == 1) 0.512 0.498 0.53

As stated on page 71 of the documentation, the use of the mean option (shortened in the code below to me), reproduces the results given by Stata’s **svy: prop** command.

svyciprop(~I(female==1), nhc, method="lo", df=degf(nhc))2.5% 97.5% I(female == 1) 0.512 0.498 0.53

You can also get the proportions of 0s, as shown below.

svyciprop(~I(female==0), nhc, method="mean")2.5% 97.5% I(female == 0) 0.488 0.474 0.5

Below is another way to get the proportions of 1s. The point of these examples is that you can use this type of syntax to get the proportion of any level of a categorical variable.

svyciprop(~I(female==1), nhc, method="mean")2.5% 97.5% I(female == 1) 0.512 0.498 0.53

Finally, let’s see how to get totals. We will use the **svytotal** function.

svytotal(~dmdborn4,design = nhc, na = TRUE)total SE dmdborn4 260493506 19670647

In the next example, we will get the coefficient of variation for this total.

cv(svytotal(~dmdborn4,design = nhc, na = TRUE))dmdborn4 dmdborn4 0.075513

Let’s get the design effect for the total.

svytotal(~dmdborn4,design = nhc, na = TRUE, deff = TRUE)total SE DEff dmdborn4 260493506 19670647 314.9

### Descriptive statistics for categorical variables

Let’s see some ways to get descriptive statistics for categorical variables, whether or not these variables are binary. We will start with the **table** function. In the example below, we get the frequencies for the variable **female**.

svytable(~female, design = nhc)female 0 1 149630839 156959842

We can also use the **table** function to get crosstabulations. We will start with two-way crosstabs.

# 2-waysvytable(~female+dmdborn4, nhc)dmdborn4 female 0 1 0 22449131 127102676 1 23543299 133390830

In the next example, we use a different syntax to do the same thing. Notice that the output is displayed differently, although the information in the output is the same.

# 2-waysvytable(~interaction(female, dmdborn4), design = nhc)interaction(female, dmdborn4) 0.0 1.0 0.1 1.1 22449131 23543299 127102676 133390830

Now let’s get a three-way table.

# 3-waysvytable(~interaction(female, dmdborn4, hsq571), design = nhc)interaction(female, dmdborn4, hsq571) 0.0.0 1.0.0 0.1.0 1.1.0 0.0.1 1.0.1 0.1.1 1.1.1 16994772.6 16101456.8 78140097.2 83984812.9 503957.9 277883.3 7083759.3 6155573.1

Let’s go a little crazy and get a four-way table.

# 4-waysvytable(~interaction(female, dmdborn4, hsq571, paq665), design = nhc)interaction(female, dmdborn4, hsq571, paq665) 0.0.0.0 1.0.0.0 0.1.0.0 1.1.0.0 0.0.1.0 1.0.1.0 0.1.1.0 9742281.10 9259759.05 41213690.37 43983640.78 314995.28 95874.61 2938039.53 1.1.1.0 0.0.0.1 1.0.0.1 0.1.0.1 1.1.0.1 0.0.1.1 1.0.1.1 2587794.24 7252491.54 6841697.75 36908800.64 40001172.12 188962.61 182008.74 0.1.1.1 1.1.1.1 4145719.79 3567778.91

Although not a descriptive statistic, let’s see how to get a chi-squared test while we are talking about tables. Of course, only a two-way table can be specified.

svychisq(~female+dmdborn4, nhc, statistic="adjWald")Design-based Wald test of association data: svychisq(~female + dmdborn4, nhc, statistic = "adjWald") F = 0.00019185, ndf = 1, ddf = 17, p-value = 0.9891

### Graphing of continuous variables

Let’s start with a histogram. By default, the density is shown on the y-axis.

svyhist(~pad630, nhc)

Instead of the density on the y-axis, we can request the count. Notice the large count of respondents in the last column on the right. The coding of the **ridageyr** variable explains this.

svyhist(~ridageyr, nhc, probability = FALSE)

We can also create boxplots.

svyboxplot(~hsq496~1, nhc, all.outliers=TRUE)

We can break the boxplot by a grouping variable. The grouping variable must be a factor.

svyboxplot(~hsq496~factor(female), nhc, all.outliers=TRUE)

We can make barcharts. In the example below, we also preview some syntax used for subpopulation analysis.

barplt<-svyby(~pad675+pad630, ~female, nhc, na = TRUE, svymean) barplot(barplt,beside=TRUE,legend=TRUE)

dotchart(barplt)

We can make a scatterplot with the sampling weights corresponding to the bubble size.

svyplot(~pad675+pad630, nhc, style="bubble")

There are a variety of density and smoothed plots that can be made. Some examples are below.

smth<-svysmooth(~pad630, design=nhc)plot(smth)

dens<-svysmooth(~pad630, design=nhc,bandwidth=30)plot(dens)

dens1<-svysmooth(~pad630, design=nhc)plot(dens1)

### Subpopulation Analysis

Before we continue with our descriptive statistics, we should pause to discuss the analysis of subpopulations. The analysis of subpopulations is one place where survey data and experimental data are quite different. If you have data from an experiment (or quasi-experiment), and you want to analyze the responses from, say, just the women, or just people over age 50, you can just delete the unwanted cases from the data set or use the by: prefix. Complex survey data are different. With survey data, you (almost) never get to delete any cases from the data set, even if you will never use them in any of your analyses. Instead, the survey package has two options that allow you to correctly analyze subpopulations of your survey data. These options are svyby and subset.survey.design. The subset.survey.design option is sort of like deleting unwanted cases (without really deleting them, of course), and the svyby option is very similar to by-group processing in that the results are shown for each group of the by-variable.

First, however, let’s take a second to see why deleting cases from a survey data set can be so problematic. There are two formulas that can used to calculate the standard errors. One formula is used when you do by-group processing or delete unwanted cases from the dataset, and survey statisticians call this the conditional approach. This is used when members of the subpopulation cannot appear in certain strata and therefore those strata should not be used in the calculation of the standard error. In practice, this rarely happens in public-use complex survey datasets. One reason is because the analyst usually does not know which combination of variables defines a particular stratum.

The other formula is used when you use the svyby option, and survey statisticians call this the unconditional approach. This is used when members of the subpopulation can be in any of the strata, even if there are some strata in the sample data that do not contain any members of the subpopulation. Because members of the subpopulation, all of the strata need to be used in the calculation of the standard error, and hence all of the data must be in the dataset. If the data set is subset (meaning that observations not to be included in the subpopulation are deleted from the data set), the standard errors of the estimates cannot be calculated correctly. When the svyby option is used, only the cases defined by the subpopulation are used in the calculation of the estimate, but all cases are used in the calculation of the standard errors. For more information on this issue, please see Sampling Techniques, Third Edition by William G. Cochran (1977) and Small Area Estimation by J. N. K. Rao (2003). A nice description of this issue given in Brady West’s 2009 Stata Conference (in Washington, D.C.).

Both **svyby** and **subset.svy.design** use the formula for the unconditional standard errors.

Let’s start by calculating the mean of **age**.

svymean(~ridageyr, nhc)mean SE ridageyr 37.185 0.6965

Now let’s calculate the mean of **age** for males and females. In this example, the variable **female** is the subpopulation variable.

svyby(~ridageyr, ~female, nhc, svymean)female ridageyr se 0 0 36.22918 0.8431945 1 1 38.09657 0.6713502

You can use more than one categorical variable to define the subpopulation. To do so, put + between the variables.

svyby(~ridageyr, ~dmdmartl+female, nhc, svymean)dmdmartl female ridageyr se 1.0 1 0 51.34608 0.7473752 2.0 2 0 67.49342 2.6062013 3.0 3 0 52.42614 0.7270806 4.0 4 0 47.16402 1.2085674 5.0 5 0 33.78962 1.3163405 6.0 6 0 40.55277 1.8945413 1.1 1 1 49.33443 0.6176356 2.1 2 1 71.39812 0.8627829 3.1 3 1 53.22067 0.7138037 4.1 4 1 47.16444 1.2200152 5.1 5 1 33.64747 1.4922117 6.1 6 1 35.94412 1.4423662

In the next example, three variables are used.

svyby(~pad630, ~dmdmartl+dmdeduc2+female, nhc, na = TRUE, svymean)dmdmartl dmdeduc2 female pad630 se 1.1.0 1 1 0 175.26250 28.12358 2.1.0 2 1 0 187.34360 62.96688 3.1.0 3 1 0 0.00000 0.00000 4.1.0 4 1 0 205.13325 35.53482 5.1.0 5 1 0 221.33212 60.76421 6.1.0 6 1 0 188.91514 71.24062 1.2.0 1 2 0 191.89233 18.18910 2.2.0 2 2 0 167.20401 14.10054 3.2.0 3 2 0 238.10555 29.71073 4.2.0 4 2 0 199.08510 38.03347 5.2.0 5 2 0 260.45360 41.06904 6.2.0 6 2 0 257.39812 43.61217 1.3.0 1 3 0 165.88490 9.75834 2.3.0 2 3 0 103.25693 20.07422 3.3.0 3 3 0 159.62009 17.52999 4.3.0 4 3 0 279.83812 48.20102 5.3.0 5 3 0 169.78003 18.53166 6.3.0 6 3 0 215.15772 24.03397 1.4.0 1 4 0 155.18399 13.21638 2.4.0 2 4 0 125.93947 35.75450 3.4.0 3 4 0 271.14720 56.91977 4.4.0 4 4 0 81.12185 29.50234 5.4.0 5 4 0 125.94122 11.38165 6.4.0 6 4 0 198.18129 25.91694 1.5.0 1 5 0 133.26699 10.84879 2.5.0 2 5 0 123.88774 52.93405 3.5.0 3 5 0 181.38341 33.00260 4.5.0 4 5 0 237.32606 84.81068 5.5.0 5 5 0 104.57660 26.15903 6.5.0 6 5 0 105.54748 7.58048 1.1.1 1 1 1 123.18726 33.99413 2.1.1 2 1 1 176.09593 72.50964 3.1.1 3 1 1 54.28887 15.49709 4.1.1 4 1 1 169.35988 80.57609 5.1.1 5 1 1 309.99124 93.15230 6.1.1 6 1 1 120.00000 0.00000 1.2.1 1 2 1 143.36162 14.00536 2.2.1 2 2 1 90.60535 17.21293 3.2.1 3 2 1 163.20222 47.28731 4.2.1 4 2 1 165.84824 46.57393 5.2.1 5 2 1 166.33346 46.45898 6.2.1 6 2 1 144.99504 24.39718 1.3.1 1 3 1 111.80803 12.92772 2.3.1 2 3 1 110.48051 37.21189 3.3.1 3 3 1 145.50738 25.43471 4.3.1 4 3 1 302.02249 117.29823 5.3.1 5 3 1 190.28528 55.64751 6.3.1 6 3 1 98.81543 21.16077 1.4.1 1 4 1 146.28879 22.89246 2.4.1 2 4 1 101.83852 11.74256 3.4.1 3 4 1 129.58833 28.49053 4.4.1 4 4 1 142.53219 26.01105 5.4.1 5 4 1 109.54923 15.62974 6.4.1 6 4 1 152.59225 23.74680 1.5.1 1 5 1 88.91933 8.30694 2.5.1 2 5 1 260.99621 63.18704 3.5.1 3 5 1 76.00663 11.32266 4.5.1 4 5 1 117.92163 34.79595 5.5.1 5 5 1 126.51918 19.29269 6.5.1 6 5 1 89.01130 23.17592

Sometimes you don’t want so much output. Rather, you just want the output for a specific group. You can get this by creating a subpopulation of the data with the subset function. In the example below, we obtain the output only for males.

smale <- subset(nhc,female == 0)summary(smale)Stratified 1 - level Cluster Sampling design (with replacement) With (31) clusters. subset(nhc, female == 0) Probabilities: Min. 1st Qu. Median Mean 3rd Qu. Max. 5.095e-06 3.044e-05 5.680e-05 6.521e-05 9.250e-05 2.326e-04 Stratum Sizes: 90 91 92 93 94 95 96 97 98 99 100 101 102 103 obs 426 496 452 288 351 355 353 289 352 338 355 357 315 129 design.PSU 3 3 3 2 2 2 2 2 2 2 2 2 2 2 actual.PSU 3 3 3 2 2 2 2 2 2 2 2 2 2 2 Data variables: [1] "dmdborn4" "dmdeduc2" "dmdmartl" "dmdyrsus" "female" "hsd010" "hsq496" [8] "hsq571" "pad615" "pad630" "pad660" "pad675" "pad680" "paq610" [15] "paq625" "paq640" "paq655" "paq665" "paq670" "paq710" "paq715" [22] "ridageyr" "sdmvpsu" "sdmvstra" "seqn" "wtint2yr"svymean(~ridageyr,design=smale)mean SE ridageyr 36.229 0.8432

### Models

A wide variety of statistical models can be run with complex survey data. With only a few exceptions, the results of these analyses can be interpreted just as the results from the same analyses with experimental or quasi-experimental data. For example, if you run an OLS regression with weighted data, assuming that the sampling plan has been correctly specified, the regression coefficients are interpreted exactly as any other OLS regression coefficient. The same is true for the various logistic regression models, including binary logistic regression, ordinal logistic regression and multinomial logistic regression (of which there is not an example in this workshop). Most of the assumptions of these models are also the same. However, some assumptions, such as the assumption regarding the normality of the residuals in OLS regression, are often not meaningful because of the large sample size commonly seen with complex survey data.

### t-tests

Let’s start doing some simple statistics. One-sample, paired-sample and independent-samples t-tests can be run with the **svyttest** function.

A one-sample t-test is shown below.

svyttest(pad675~0, nhc, na = TRUE)Design-based one-sample t-test data: pad675 ~ 0 t = 40.447, df = 16, p-value < 2.2e-16 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 63.84046 70.90254 sample estimates: mean 67.3715

In the example of the paired-samples t-test below, the “I” is used to tell R to leave the part in parentheses “as is”, meaning do the subtraction between the two variables. Hence, the formula means that the difference between

svyttest(I(pad660-pad675)~0, nhc, na = TRUE)Design-based one-sample t-test data: I(pad660 - pad675) ~ 0 t = 3.3059, df = 16, p-value = 0.004464 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 2.718683 12.437915 sample estimates: mean 7.578299

Now let’s run an independent-samples t-test.

svyttest(pad630~female, nhc)Design-based t-test data: pad630 ~ female t = -5.6277, df = 16, p-value = 3.779e-05 alternative hypothesis: true difference in mean is not equal to 0 95 percent confidence interval: -45.76433 -22.12149 sample estimates: difference in mean -33.94291

As you probably know, an independent-samples t-test tests the null hypothesis that the difference in the means of the two groups is 0. Another way to think about this type of t-test is to think of it as a linear regression with a single binary predictor. The intercept will be the mean of the reference group, and the coefficient will be the difference between the two groups.

We will start by running the **t-tes**t function as before, and then replicate the results using the **svyglm function**, which can be used to run a linear regression. The **svyby** function is used with the **covmat** argument to save the elements to a matrix so that we can use the **svycontrast** function to subtract the values. The purpose of this example is not to belabor the point about a t-test, but rather to show how to get a matrix of values and then compare those values with the **svycontrast** function in a simple example where the answer is already known.

svyttest(ridageyr~female, nhc)Design-based t-test data: ridageyr ~ female t = 2.9691, df = 16, p-value = 0.009043 alternative hypothesis: true difference in mean is not equal to 0 95 percent confidence interval: 0.634692 3.100076 sample estimates: difference in mean 1.867384

summary(svyglm(ridageyr~female, design=nhc))Call: svyglm(formula = ridageyr ~ female, design = nhc) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 36.2292 0.8432 42.967 < 2e-16 *** female 1.8674 0.6289 2.969 0.00904 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 499.5325) Number of Fisher Scoring iterations: 2

a <- svyby(~ridageyr, ~female, nhc, na.rm.by = TRUE, svymean, covmat = TRUE)vcov(a)0 1 0 0.7109769 0.3830638 1 0.3830638 0.4507111 a female ridageyr se 0 0 36.22918 0.8431945 1 1 38.09657 0.6713502

svycontrast(a, c( -1, 1))contrast SE contrast 1.8674 0.6289 # 4.589723 - 6.153479

This example is similar to the previous example, except that here the **svypredmeans** function is used.

summary(svyglm(pad630~female, design=nhc))Call: svyglm(formula = pad630 ~ female, design = nhc) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 155.627 7.008 22.206 1.90e-13 *** female -33.943 6.031 -5.628 3.78e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 21951.7) Number of Fisher Scoring iterations: 2

Note that the variable **female** cannot be in model if you want to get the predicted means for that variable.

ttest1 <- (svyglm(ridageyr~1, design=nhc))summary(ttest1) Call: svyglm(formula = ridageyr ~ 1, design = nhc) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 37.1852 0.6965 53.39 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 500.4039) Number of Fisher Scoring iterations: 2

The variable **female** is used here to get the predicted means for each level of **female**.

svypredmeans(ttest1, ~female)mean SE 0 36.229 0.8432 1 38.097 0.6714

tt<-svyttest(pad630~female, nhc) ttDesign-based t-test data: pad630 ~ female t = -5.6277, df = 16, p-value = 3.779e-05 alternative hypothesis: true difference in mean is not equal to 0 95 percent confidence interval: -45.76433 -22.12149 sample estimates: difference in mean -33.94291

We can get the confidence interval around the difference. In this example, we get the 90% confidence interval.

confint(tt, level=0.9)[1] -43.67864 -24.20718 attr(,"conf.level") [1] 0.9

### Multiple linear regression

We need to use the **summary** function to get the standard errors, test statistics and p-values. Let’s start with a model that has no interaction terms. The outcome variable will be **pad630**, and the predictors will be **female** and **hsq571**.

summary(svyglm(pad630~female+hsq571, design=nhc, na.action = na.omit))Call: svyglm(formula = pad630 ~ female + hsq571, design = nhc, na.action = na.omit) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 164.707 6.595 24.975 1.24e-13 *** female -39.296 6.325 -6.212 1.66e-05 *** hsq571 -11.722 20.267 -0.578 0.572 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 23558.77) Number of Fisher Scoring iterations: 2

Now let’s add an interaction between the two predictor variables.

summary(svyglm(pad630~female*hsq571, design=nhc, na.action = na.omit))Call: svyglm(formula = pad630 ~ female * hsq571, design = nhc, na.action = na.omit) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 167.107 6.847 24.404 7.13e-13 *** female -44.520 6.827 -6.521 1.35e-05 *** hsq571 -40.975 18.282 -2.241 0.04174 * female:hsq571 67.630 19.484 3.471 0.00374 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 23441.51) Number of Fisher Scoring iterations: 2

glm1 <- (svyglm(pad630~female*hsq571, design=nhc, na.action = na.omit))glm1Stratified 1 - level Cluster Sampling design (with replacement) With (31) clusters. svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Call: svyglm(formula = pad630 ~ female * hsq571, design = nhc, na.action = na.omit) Coefficients: (Intercept) female hsq571 female:hsq571 167.11 -44.52 -40.97 67.63 Degrees of Freedom: 1672 Total (i.e. Null); 14 Residual (8083 observations deleted due to missingness) Null Deviance: 40340000 Residual Deviance: 39190000 AIC: 21580 confint(glm1) 2.5 % 97.5 % (Intercept) 153.68606 180.527513 female -57.89966 -31.139511 hsq571 -76.80715 -5.141942 female:hsq571 29.44252 105.817225

This example is just like the previous one, only here factor notation is used. This is important when the categorical predictor has more than two levels.

summary(svyglm(pad630~factor(female)*factor(dmdmartl), design=nhc, na.action = na.omit))Call: svyglm(formula = pad630 ~ factor(female) * factor(dmdmartl), design = nhc, na.action = na.omit) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 156.694 8.766 17.876 1.97e-06 *** factor(female)1 -38.069 12.438 -3.061 0.0222 * factor(dmdmartl)2 -24.083 23.185 -1.039 0.3390 factor(dmdmartl)3 52.729 28.804 1.831 0.1169 factor(dmdmartl)4 31.274 32.355 0.967 0.3711 factor(dmdmartl)5 -10.271 13.330 -0.771 0.4702 factor(dmdmartl)6 44.728 16.959 2.637 0.0387 * factor(female)1:factor(dmdmartl)2 26.414 28.913 0.914 0.3962 factor(female)1:factor(dmdmartl)3 -46.720 43.539 -1.073 0.3245 factor(female)1:factor(dmdmartl)4 51.033 51.334 0.994 0.3585 factor(female)1:factor(dmdmartl)5 23.355 18.086 1.291 0.2441 factor(female)1:factor(dmdmartl)6 -40.780 25.360 -1.608 0.1589 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 23215.65) Number of Fisher Scoring iterations: 2ols1 <- (svyglm(pad630~1, design=nhc, na.action = na.omit)) predmarg<-svypredmeans(ols1, ~interaction(female,dmdmartl)) predmargmean SE 1.1 118.63 9.5904 0.6 201.42 13.4189 0.1 156.69 8.7656 0.5 146.42 12.3533 1.5 131.71 13.3056 1.4 200.93 52.4999 1.6 122.57 14.8851 0.3 209.42 27.2862 1.2 120.96 11.1044 1.3 124.63 15.5059 0.4 187.97 38.8611 0.2 132.61 20.7758

### Non-parametric tests

Non-parametric tests can also be done. Let’s start with a Wilcoxon signed rank test, which is the non-parametric analog of an independent-samples t-test.

wil <- svyranktest(hsq496~female, design = nhc, na = TRUE, test = c("wilcoxon"))wilDesign-based KruskalWallis test data: hsq496 ~ female t = 6.3291, df = 16, p-value = 1.002e-05 alternative hypothesis: true difference in mean rank score is not equal to 0 sample estimates: difference in mean rank score 0.06896535

This is an example of a median test.

mtest <- svyranktest(hsq496~female, design = nhc, na = TRUE, test=("median"))mtestDesign-based median test data: hsq496 ~ female t = 4.9504, df = 16, p-value = 0.0001446 alternative hypothesis: true difference in mean rank score is not equal to 0 sample estimates: difference in mean rank score 0.11726

This is an example of a Kruskal Wallis test, which is the non-parametric analog of a one-way ANOVA.

kwtest <- svyranktest(hsq496~female, design = nhc, na = TRUE, test=("KruskalWallis"))kwtestDesign-based KruskalWallis test data: hsq496 ~ female t = 6.3291, df = 16, p-value = 1.002e-05 alternative hypothesis: true difference in mean rank score is not equal to 0 sample estimates: difference in mean rank score 0.06896535

### Logistic regression

Let’s see a few examples of logistic regression.

logit1 <- (svyglm(paq665~factor(hsd010)+ridageyr, family=quasibinomial, design=nhc, na.action = na.omit))summary(logit1)Call: svyglm(formula = paq665 ~ factor(hsd010) + ridageyr, design = nhc, family = quasibinomial, na.action = na.omit) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.717830 0.123638 5.806 8.40e-05 *** factor(hsd010)2 -0.053653 0.099863 -0.537 0.600902 factor(hsd010)3 -0.541820 0.104759 -5.172 0.000232 *** factor(hsd010)4 -0.981956 0.103744 -9.465 6.46e-07 *** factor(hsd010)5 -1.882124 0.201158 -9.356 7.31e-07 *** ridageyr -0.009443 0.002191 -4.311 0.001013 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasibinomial family taken to be 1.219511) Number of Fisher Scoring iterations: 4

In the next example, we will run the logistic regression on a subpopulation (respondents over age 20).

subset1 <- subset(nhc, ridageyr > 20)logit2 <- (svyglm(paq665~factor(hsd010)+ridageyr, family=quasibinomial, design=subset1, na.action = na.omit))summary(logit2) Call: svyglm(formula = paq665 ~ factor(hsd010) + ridageyr, design = subset1, family = quasibinomial, na.action = na.omit) Survey design: subset(nhc, ridageyr > 20) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.86401 0.17432 4.956 0.000333 *** factor(hsd010)2 -0.10723 0.11507 -0.932 0.369739 factor(hsd010)3 -0.59868 0.11209 -5.341 0.000176 *** factor(hsd010)4 -1.09325 0.11226 -9.738 4.77e-07 *** factor(hsd010)5 -2.04923 0.19386 -10.571 1.96e-07 *** ridageyr -0.01097 0.00292 -3.756 0.002740 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasibinomial family taken to be 1.024088) Number of Fisher Scoring iterations: 4

We can also get a Wald test for a variable in the model.

regTermTest(logit2, ~ridageyr)Wald test for ridageyr in svyglm(formula = paq665 ~ factor(hsd010) + ridageyr, design = subset1, family = quasibinomial, na.action = na.omit) F = 14.10816 on 1 and 12 df: p= 0.0027404

Instead of getting an R-squared value as you do in linear regression, a pseudo-R-squared is given in logistic regression. There are many different versions of pseudo-R-squared, and two of them are available with the **psrsq** function.

psrsq(logit2, method = c("Cox-Snell"))[1] 0.05148869psrsq(logit2, method = c("Nagelkerke"))[1] 0.06873682

### Ordered logistic regression

Below is an example of an ordered logistic regression. Note that the outcome variable must be a factor.

ologit1 <- svyolr(factor(dmdeduc2)~factor(female)+factor(dmdborn4)+pad680, design = nhc, method = c("logistic"))summary(ologit1)Call: svyolr(factor(dmdeduc2) ~ factor(female) + factor(dmdborn4) + pad680, design = nhc, method = c("logistic")) Coefficients: Value Std. Error t value factor(female)1 0.097569476 0.0487417745 2.001763 factor(dmdborn4)1 0.709389138 0.1215114128 5.838045 pad680 0.001923994 0.0001494462 12.874160 Intercepts: Value Std. Error t value 1|2 -1.5280 0.1541 -9.9183 2|3 -0.3101 0.1562 -1.9849 3|4 0.8095 0.1344 6.0214 4|5 2.2157 0.1210 18.3118 (4234 observations deleted due to missingness)

## Poisson regression

Poisson regression can be run. This is a type of count model (meaning that the outcome variable should be a count).

summary(svyglm(pad675~female, design=nhc, family=poisson()))Call: svyglm(formula = pad675 ~ female, design = nhc, family = poisson()) Survey design: svydesign(id = ~sdmvpsu, weights = ~wtint2yr, strata = ~sdmvstra, nest = TRUE, survey.lonely.psu = "adjust", data = nhanes2012) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.34112 0.03487 124.491 < 2e-16 *** female -0.27148 0.03698 -7.341 1.66e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 66.18153) Number of Fisher Scoring iterations: 5

### Other types of analyses available in the survey package

There are many more types of analyses that are available in the **survey** package and in other packages that work with complex survey data. Below are a few examples.

The example below shows a principle components analysis (PCA).

pc <- svyprcomp(~pad630+pad675+hsd010, design=nhc,scale=TRUE,scores=TRUE)pcStandard deviations (1, .., p=3): [1] 1.3023573 0.8183481 0.7963491 Rotation (n x k) = (3 x 3): PC1 PC2 PC3 pad630 0.5769018 -0.5906013 -0.5642468 pad675 0.5690424 0.7861726 -0.2410879 hsd010 0.5859822 -0.1819963 0.7896216biplot(pc, weight="scaled")

This is an example of Cronbach’s alpha.

svycralpha(~hsq571+dmdborn4, design=nhc, na.rm = TRUE)*alpha* 0.1339271