As a statistical programming language, R allows users to access precise statistics instead of simply printing a mass of output to the screen. The examples below highlight how to create a complex sample survey design object and then directly query specific coefficients, error terms, and other survey design-related information as needed.

This page uses the `survey`

package. Make sure that you can load
it before trying to run the examples on this page. If you do not have
a package installed, run: `install.packages("packagename")`

, or
if you see the version is out of date, run: `update.packages()`

.

```
library(survey)
```

**Version info: **Code for this page was tested in R version 3.0.1 (2013-05-16)
On: 2013-06-25
With: survey 3.29-5; knitr 1.2

## Example

This example is taken from Lehtonen and Pahkinen’s Practical Methods for Design and Analysis of Complex Surveys. Page 46 Table 2.6 Estimates from a systematic sample drawn from the Province’91 population using implicit stratification.

NOTE: The standard error of the total is different from that shown in the text (the text shows 11802). However, we get the 13627 in each of the statistical packages in which we have tried to recreate this example.

Import the dataset from text directly into R using the `read.table`

function and the `text=`

parameter specifying the entire data set. The syntax `n`

indicates the end of one line of data.

```
mydata <-
read.table( text =
"id str clu wt ue91 lab91
1 1 1 4 4123 33786
2 1 5 4 721 4930
3 2 9 4 194 2069
4 2 13 4 129 927
5 2 17 4 239 2144
6 2 21 4 61 573
7 2 25 4 262 1737
8 2 29 4 166 1615" ,
header = TRUE
)
```

More detail about `read.table`

or any other R function can be viewed by typing a question mark in front of the function name in the R console. For example: `?read.table`

.

In the R language, individual data sets are stored as `data.frame`

objects, allowing users to load as many tables into working memory as necessary for the analysis. After loading the `mydata`

table into memory, R functions can be run directly on this data table.

# the `class` of the `mydata` object class(mydata)

## [1] "data.frame"

# the first six rows of the data head(mydata)

## id str clu wt ue91 lab91 ## 1 1 1 1 4 4123 33786 ## 2 2 1 5 4 721 4930 ## 3 3 2 9 4 194 2069 ## 4 4 2 13 4 129 927 ## 5 5 2 17 4 239 2144 ## 6 6 2 21 4 61 573

# the last six rows of the data tail(mydata)

## id str clu wt ue91 lab91 ## 3 3 2 9 4 194 2069 ## 4 4 2 13 4 129 927 ## 5 5 2 17 4 239 2144 ## 6 6 2 21 4 61 573 ## 7 7 2 25 4 262 1737 ## 8 8 2 29 4 166 1615

# the number of columns ncol(mydata)

## [1] 6

View a simple summary of your variable of interest. This `summary`

function accesses the `ue91`

column (variable) stored inside the `mydata`

data.frame object.

```
summary(mydata$ue91)
```

## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 61 157 216 737 377 4120

Initiate your `svydesign`

object for a **systematic sample** design. This ``mydesign``

object will be used for all subsequent analysis commands:

```
mydesign <-
svydesign(
id = ~clu ,
data = mydata ,
weight = ~wt ,
strata = ~str
)
```

From this point forward, the sampling specifications of the `mydata`

data set’s survey design have been fixed and most analysis commands will simply use the set of tools outlined on the R survey package homepage, referring to the object ``mydesign``

at the `design=`

parameter of the specific R function or method.

View the survey design structure:

mydesign

## Stratified Independent Sampling design (with replacement) ## svydesign(id = ~clu, data = mydata, weight = ~wt, strata = ~str)

View the survey design’s object class or type:

```
class(mydesign)
```

## [1] "survey.design2" "survey.design"

View the weighted total population of this survey design, by referring to the `wt`

column from the original data.frame:

```
sum(mydata$wt)
```

## [1] 32

View the number of unique PSUs (clusters) in this survey design, by referring to the `clu`

column from the original data.frame:

length(unique(mydata$clu))

## [1] 8

View the number of unique strata in this survey design, by referring to the `str`

column from the original data.frame:

length(unique(mydata$str))

## [1] 2

View the degrees of freedom for this survey design object:

```
degf(mydesign)
```

## [1] 6

Count the number of unweighted observations where the variable `ue91`

is not missing:

```
unwtd.count(~ue91, mydesign)
```

## counts SE ## counts 8 0

Print the mean and standard error of the `ue91`

variable:

```
svymean(~ue91, mydesign)
```

## mean SE ## ue91 737 426

Print the weighted total and standard error of the `ue91`

variable:

```
svytotal(~ue91, mydesign)
```

## total SE ## ue91 23580 13627

Alternatively, the result of a function call like `svymean`

or `svytotal`

can be stored into a secondary R object.

mysvymean <- svymean(~ue91, mydesign, deff = TRUE)

Once created, this `svymean`

can be queried independently from the `mydesign`

object. For example, the `coef`

and `SE`

functions can directly extract those attributes:

```
coef(mysvymean)
```

## ue91 ## 737

```
SE(mysvymean)
```

## ue91 ## ue91 426

The design effect extraction function `deff`

can only be used if the original `svymean`

call that created the object `mysvymean`

included the parameter `deff = TRUE`

.

```
deff(mysvymean)
```

## ue91 ## 1.01

Also note that the number of decimal places shown can be adjusted by modifying the `digits`

parameter within the `options`

function at any time.

options(digits = 10) SE(mysvymean)

## ue91 ## ue91 425.8496397

We can use the `confint`

function to obtain confidence intervals for the coefficient estimates.

```
confint(mysvymean)
```

## 2.5 % 97.5 % ## ue91 -97.77495662 1571.524957

Save the ratio of `ue91`

to `lab91`

into a new object `myratio`

and at the same time print it to the screen by encapsulaing the entire statement in parentheses.

```
(myratio <- svyratio(~ue91, ~lab91, mydesign))
```

## Ratio estimator: svyratio.survey.design2(~ue91, ~lab91, mydesign) ## Ratios= ## lab91 ## ue91 0.1233754003 ## SEs= ## lab91 ## ue91 0.003848016044

We can then use the `confint`

function to obtain confidence intervals for the ratio.

```
confint(myratio)
```

## 2.5 % 97.5 % ## ue91/lab91 0.1158334274 0.1309173731

We can specify the `df=`

parameter to use the survey design’s degrees of freedom (instead of the default `df=Inf`

) to replicate Stata’s confidence interval calculation method.

# matches stata confint(myratio, df = degf(mydesign))

## 2.5 % 97.5 % ## ue91/lab91 0.1139596442 0.1327911563

Print the median of the `ue91`

variable, including the confidence interval in the output.

svyquantile(~ue91, mydesign, quantiles = 0.5, ci = TRUE)

## $quantiles ## 0.5 ## ue91 194 ## ## $CIs ## , , ue91 ## ## 0.5 ## (lower 96.41575121 ## upper) 481.94367931

## See also

- The R survey package homepage
- Lumley, T. Complex Surveys: A Guide to Analysis Using R (Wiley Series in Survey Methodology)
- Damico, A. Step-by-step instructions to analyze major public-use survey data sets with the R language