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 following packages. Make sure that you can load
them 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(foreign) 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; foreign 0.8-54; knitr 1.2
Example 1
This example is taken from Levy and Lemeshow’s Sampling of Populations.
Page 350 cluster sampling with unequal probabilities: probability proportional to size sampling
Import the Stata dataset directly into R using the read.dta function from the foreign package:
mydata <- read.dta( "https://stats.idre.ucla.edu/stat/examples/sop/hospslct.dta" , convert.factors = FALSE )
More detail about read.dta 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.dta.
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)
## drawing hospno admiss lifethrt dxdead p1 p2 p wstar ## 1 1 2 5036 1 1 0.411 0.00198 0.000816 1001 ## 2 1 2 5036 0 0 0.411 0.00198 0.000816 1001 ## 3 1 2 5036 0 0 0.411 0.00198 0.000816 1001 ## 4 1 2 5036 0 0 0.411 0.00198 0.000816 1001 ## 5 1 2 5036 0 0 0.411 0.00198 0.000816 1001 ## 6 1 2 5036 0 0 0.411 0.00198 0.000816 1001
# the last six rows of the data tail(mydata)
## drawing hospno admiss lifethrt dxdead p1 p2 p wstar ## 45 5 9 4672 0 0 0.387 0.00214 0.000828 1001 ## 46 5 9 4672 0 0 0.387 0.00214 0.000828 1001 ## 47 5 9 4672 0 0 0.387 0.00214 0.000828 1001 ## 48 5 9 4672 0 0 0.387 0.00214 0.000828 1001 ## 49 5 9 4672 0 0 0.387 0.00214 0.000828 1001 ## 50 5 9 4672 0 0 0.387 0.00214 0.000828 1001
# the number of columns ncol(mydata)
## [1] 9
View a simple tabulation of your variable of interest. This table function accesses the lifethrt column (variable) stored inside the mydata data.frame object.
table(mydata$lifethrt)
## ## 0 1 ## 44 6
Initiate your svydesign object for a probability proportional to size sampling design. This `mydesign` object will be used for all subsequent analysis commands:
mydesign <-
svydesign(
id = ~drawing ,
data = mydata ,
weight = ~wstar
)
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
## 1 - level Cluster Sampling design (with replacement) ## With (5) clusters. ## svydesign(id = ~drawing, data = mydata, weight = ~wstar)
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 wstar column from the original data.frame:
sum(mydata$wstar)
## [1] 50056
View the number of unique PSUs (clusters) in this survey design, by referring to the drawing column from the original data.frame:
length(unique(mydata$drawing))
## [1] 5
View the degrees of freedom for this survey design object:
degf(mydesign)
## [1] 4
Count the number of unweighted observations where the variable lifethrt is not missing:
unwtd.count(~lifethrt, mydesign)
## counts SE ## counts 50 0
Print the mean and standard error of two variables:
svymean(~lifethrt + dxdead, mydesign)
## mean SE ## lifethrt 0.12 0.02 ## dxdead 0.04 0.02
Print the weighted total and standard error of the same two variables:
svytotal(~lifethrt + dxdead, mydesign)
## total SE ## lifethrt 6007 1001 ## dxdead 2002 1226
Alternatively, the result of a function call like svymean or svytotal can be stored into a secondary R object.
mysvymean <- svymean(~lifethrt + dxdead, 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)
## lifethrt dxdead ## 0.12 0.04
SE(mysvymean)
## lifethrt dxdead ## 0.0200 0.0245
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)
## lifethrt dxdead ## 0.186 0.766
We can use the confint function to obtain confidence intervals for the coefficient estimates.
confint(mysvymean)
## 2.5 % 97.5 % ## lifethrt 0.08080 0.159 ## dxdead -0.00801 0.088
Note from our unweighted table, the variable lifethrt was binary (composed strictly of zeroes and ones). To produce confidence intervals more accurate for proportions, users might start with the options discussed in ?svyciprop. For example:
svyciprop(~lifethrt, mydesign, method = "logit")
## 2.5% 97.5% ## lifethrt 0.1200 0.0898 0.16
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)
## lifethrt dxdead ## 0.02000000000 0.02449489743
Save the ratio of dxdead to lifethrt into a new object myratio and at the same time print it to the screen by encapsulaing the entire statement in parentheses.
( myratio <- svyratio( ~dxdead , ~lifethrt , mydesign ) )
## Ratio estimator: svyratio.survey.design2(~dxdead, ~lifethrt, mydesign) ## Ratios= ## lifethrt ## dxdead 0.3333333333 ## SEs= ## lifethrt ## dxdead 0.2324055629
We can then use the confint function to obtain confidence intervals for the ratio.
confint(myratio)
## 2.5 % 97.5 % ## dxdead/lifethrt -0.1221731998 0.7888398665
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 % ## dxdead/lifethrt -0.3119279543 0.9785946209
Example 2
This example is taken from Levy and Lemeshow’s Sampling of Populations.
Page 351 cluster sampling with unequal probabilities: probability proportional to size sampling
NOTE: The following recodes are necessary to create a new weight variable, w2star.
# records where `hospno` is two, create a new variable `tl` with the value 785 throughout. mydata[ mydata$hospno == 2 , 'tl' ] <- 785 # records where `hospno` is five, set `tl` to 3404 mydata[ mydata$hospno == 5 , 'tl' ] <- 3404 # records where `hospno` is nine, set `tl` to 778 mydata[ mydata$hospno == 9 , 'tl' ] <- 778 # create a new `w2star` variable using `transform`. mydata <- transform( mydata , w2star = ( admiss / 50 ) * ( 7087 / tl ) )
The variable w2star has been set equal to the product of admiss variable divided by fifty and 7,087 divided by the newly-created tl variable.
Initiate your svydesign object for a probability proportional to size sampling design. This `mydesign2` object will be used for all subsequent analysis commands:
mydesign2 <-
svydesign(
id = ~drawing ,
data = mydata ,
weight = ~w2star
)
This new survey design object differs from the mydesign object in its weighting by the newly-created w2star variable.
Print the mean and standard error of two variables:
svymean(~lifethrt + dxdead, mydesign2)
## mean SE ## lifethrt 0.121904297 0.02142 ## dxdead 0.034287108 0.02300
Print the weighted total and standard error of the same two variables:
svytotal(~lifethrt + dxdead, mydesign2)
## total SE ## lifethrt 6259.1758 1277.3220 ## dxdead 1760.4715 1079.0434
Print the ratio of dxdead to lifethrt using the new survey design object:
svyratio(~dxdead, ~lifethrt, mydesign2)
## Ratio estimator: svyratio.survey.design2(~dxdead, ~lifethrt, mydesign2) ## Ratios= ## lifethrt ## dxdead 0.2812625046 ## SEs= ## lifethrt ## dxdead 0.2114834755
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