\documentclass{beamer}\usepackage[]{graphicx}\usepackage[]{xcolor}
% maxwidth is the original width if it is less than linewidth
% otherwise use linewidth (to make sure the graphics do not exceed the margin)
\makeatletter
\def\maxwidth{ %
  \ifdim\Gin@nat@width>\linewidth
    \linewidth
  \else
    \Gin@nat@width
  \fi
}
\makeatother

\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
\let\hlipl\hlkwb

\usepackage{framed}
\makeatletter
\newenvironment{kframe}{%
 \def\at@end@of@kframe{}%
 \ifinner\ifhmode%
  \def\at@end@of@kframe{\end{minipage}}%
  \begin{minipage}{\columnwidth}%
 \fi\fi%
 \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
 \colorbox{shadecolor}{##1}\hskip-\fboxsep
     % There is no \\@totalrightmargin, so:
     \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
 \MakeFramed {\advance\hsize-\width
   \@totalleftmargin\z@ \linewidth\hsize
   \@setminipage}}%
 {\par\unskip\endMakeFramed%
 \at@end@of@kframe}
\makeatother

\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX

\usepackage{alltt}
\usetheme{Frankfurt}
\usecolortheme{orchid}

\usepackage{amsbsy}
\usepackage{hyperref}  % hyperlink
\usepackage{multimedia}
\usepackage{beamerthemesplit}
\usepackage{tikz}
\usepackage[utf8]{inputenc} 
\usepackage{amsmath,amssymb}
\usepackage{colortbl}
\usepackage[english]{babel}
\usepackage{tabularx}
\usepackage{pgf}  
\usepackage{textpos}
\usepackage{multicol}
\usepackage{chngcntr}
\usepackage{lscape}
\usepackage{amsfonts}
\usepackage{graphicx}  
\usepackage{color}
\usepackage{textcomp}
%\usepackage{wasysym}
\usepackage{float}
\usepackage{epstopdf}
\usepackage{pdfpages}


\newcommand{\E}{\mbox{E}\,}
\newcommand{\Cov}{\mbox{Cov}\,}
\newcommand{\Var}{\mbox{Var}\,}
\newcommand{\MSE}{\mbox{MSE}\,}
\newcommand{\tr}{\mbox{tr}\,}
\newcommand{\diag}{\mbox{diag}\,}
\newcommand{\bb}{{\mathbf b}}
\newcommand{\bc}{{\mathbf c}}
\newcommand{\bd}{{\mathbf d}}
\newcommand{\be}{{\mathbf e}}
\newcommand{\bff}{{\mathbf f}}
\newcommand{\bh}{{\mathbf h}}
\newcommand{\bm}{{\mathbf m}}
\newcommand{\bn}{{\mathbf n}}
\newcommand{\bt}{{\mathbf t}}
\newcommand{\bu}{{\mathbf u}}
\newcommand{\bv}{{\mathbf v}}
\newcommand{\bx}{{\mathbf x}}
\newcommand{\bw}{{\mathbf w}}
\newcommand{\bg}{{\mathbf g}}
\newcommand{\bp}{{\mathbf p}}
\newcommand{\by}{{\mathbf y}}
\newcommand{\bz}{{\mathbf z}}
\newcommand{\bA}{{\mathbf A}}
\newcommand{\bB}{{\mathbf B}}
\newcommand{\bC}{{\mathbf C}}
\newcommand{\bD}{{\mathbf D}}
\newcommand{\bE}{{\mathbf E}}
\newcommand{\bF}{{\mathbf F}}
\newcommand{\bG}{{\mathbf G}}
\newcommand{\bH}{{\mathbf H}}
\newcommand{\bI}{{\mathbf I}}
\newcommand{\bJ}{{\mathbf J}}
\newcommand{\bK}{{\mathbf K}}
\newcommand{\bL}{{\mathbf L}}
\newcommand{\bM}{{\mathbf M}}
\newcommand{\bN}{{\mathbf N}}
\newcommand{\bP}{{\mathbf P}}
\newcommand{\bQ}{{\mathbf Q}}
\newcommand{\bR}{{\mathbf R}}
\newcommand{\bS}{{\mathbf S}}
\newcommand{\bT}{{\mathbf T}}
\newcommand{\bU}{{\mathbf U}}
\newcommand{\bV}{{\mathbf V}}
\newcommand{\bW}{{\mathbf W}}
\newcommand{\bX}{{\mathbf X}}
\newcommand{\bY}{{\mathbf Y}}
\newcommand{\bZ}{{\mathbf Z}}
\newcommand{\bone}{{\mathbf 1}}
\newcommand{\bzero}{{\mathbf 0}}
\newcommand{\bal}{{\boldsymbol \alpha}}
\newcommand{\beps}{{\boldsymbol \epsilon}}
\newcommand{\vj}{v_{\tiny\mbox{JACK}}}
\newcommand{\vaj}{v_{\tiny\mbox{JACK}}^{\tiny\mbox{adj}}}
\newcommand{\stat}{\overline{\mbox{QOL4-QOL0}}}

\newcommand{\bbe}{{\boldsymbol \beta}}
\newcommand{\bga}{{\boldsymbol \gamma}}
\newcommand{\bsig}{{\boldsymbol \sigma}}
\newcommand{\bde}{{\boldsymbol \delta}}
\newcommand{\bphi}{{\boldsymbol \phi}}
\newcommand{\bpsi}{{\boldsymbol \psi}}
\newcommand{\btau}{{\boldsymbol \tau}}
\newcommand{\bDe}{{\boldsymbol \Delta}}
\newcommand{\boldeta}{{\boldsymbol \eta}}
\newcommand{\bpi}{{\boldsymbol p}}
\newcommand{\boldnabla}{{\boldsymbol \nabla}}
\newcommand{\blam}{{\boldsymbol \lambda}}
\newcommand{\bmu}{{\boldsymbol \mu}}
\newcommand{\bth}{{\boldsymbol \theta}}
\newcommand{\bthhat}{{\bf \hat {\boldsymbol \theta}}}
\newcommand{\bGam}{{\mathbf \Gamma}}
\newcommand{\bLam}{{\mathbf \Lambda}}
\newcommand{\bSig}{{\mathbf \Sigma}}
\newcommand{\bOm}{{\mathbf \Omega}}
\newcommand{\bPsi}{{\mathbf \Psi}}
\newcommand{\bUp}{{\mathbf \Upsilon}}
\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\ga}{\gamma}
\newcommand{\la}{\lambda}
\newcommand{\eps}{\epsilon}
\renewcommand{\th}{\theta}
\newcommand{\om}{\omega}
\newcommand\cF{\mathcal{F}}
\newcommand\cL{\mathcal{L}}
\newcommand\cS{\mathcal{S}}
\newcommand\cU{\mathcal{U}}
\newcommand\ep{{\varepsilon}}
\newcommand\s{{\sigma}}
\newcommand{\yhat}{\hat{Y}}
\newcommand\uth{{\textrm{th}}}
\newcommand{\ds}{\displaystyle}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}


\setbeamertemplate{footline}[frame number]
\setlength{\emergencystretch}{2pt}

\pgfdeclareimage[height=1.8cm]{logoUNM}{./plots/unmlog}
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}
\setbeamertemplate{itemize items}[triangle]
\setbeamertemplate{itemize subitem}[square]
\setbeamertemplate{itemize subsubitem}[circle]
\setbeamertemplate{enumerate items}[default]
\setbeamertemplate{enumerate subitem}[square]
\setbeamertemplate{enumerate subsubitem}[circle]
\setbeamertemplate{caption}[numbered]


\title{Stat472/572 Sampling: Design and Analysis}
\author {Instructor: Yan Lu}
\date{}

% very important to use option [fragile] for frames containing code output!

\begin{frame}
 \titlepage
 \vspace{-1cm}
\centering
    \textbf{\Large Chapter 7: Complex Surveys}
 \begin{textblock*}{2.0cm}(0.0cm,-3.0cm)
   %   \rule{2cm}{1cm}
   \pgfuseimage{logoUNM}
 \end{textblock*}
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Topic Overviews}
 \bi
 \item Complex surveys
 \item Survey plots
 \item Sampling and experimental
design
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Building blocks for surveys}

\bi
\item Cluster sampling with replacement
\item Cluster sampling without replacement
\item Stratification
\item Ratio estimation
\item Weights
\item Computer intensive methods
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Cluster sampling with replacement}
\bi
 \item Select a sample of $n$
clusters with replacement\\
----suppose cluster $i$ is selected with
probability $\psi_i$
\item Estimate the total for cluster $i$ by an
unbiased estimate $\hat t_i$
\item Estimator for population total is $\hat t_{\psi} = \ds\frac 1 n \sum\limits_{i=1}^N
\sum\limits_{j=1}^{Q_i} \ds\frac{\hat t_{ij}}{\psi_i}$
\item $\hat V(\hat t_{\psi}) = \ds\frac 1 n \sum\limits_{i=1}^N
\sum\limits_{j=1}^{Q_i} \left(\ds\frac {\hat t_{ij}}{\psi_i} -\hat
t_{\psi}\right)^2/(n-1) $
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Cluster sampling without replacement}
\bi
\item Select a sample of $n$
clusters without replacement
\item Suppose cluster $i$ is selected in the
sample with probability $\pi_i$
\item Estimate the total for cluster
$i$ by an unbiased estimate $\hat t_i$
\item Use the Horvitz-Thompson estimate of the population total
{\large
$$\hat t_{\mbox{HT}} =\sum_{i\in S} \hat t_i /\pi_i$$
\item $\hat V(\hat t_{\mbox{HT}}) = \left(1-\ds\frac n N\right) \ds\frac n {n-1}
\sum\limits_{i\in S} \left(\ds\frac{\hat t_i}{\pi_i} -
 \ds\frac 1 n \sum\limits_{j\in S} \ds\frac{\hat t_j}{\pi_j}\right)^2 + \sum\limits_{i\in S}
\ds\frac{\hat V(\hat t_i)}{\pi_i}$}
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Stratification}
\bi
\item Estimate the strata totals by $\hat t_1, \hat t_2, \cdots,\hat t_H$
\item The estimated variances for the strata totals are \\$\hat V(\hat t_1),
\hat V(\hat t_2), \cdots , \hat V(\hat t_H)$
\item The estimator of
the population total is $\hat t = \sum_{i=1}^H \hat t_i$
\item The estimator
of the variance is $\hat V(\hat t) =\sum_{i=1}^H \hat V(\hat t_i)$
\ei

\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}

{\bf Ratio estimation}
\bi
 \item Let $\hat t_x$ and $\hat t_y$ be
estimators of $t_x$ and $t_y$, respectively
\item The ratio is estimated by $\hat B=\hat t_y/\hat t_x$
\item The ratio estimator of the population total is $\hat t_{yr} = \hat B t_x$
\item The estimated variance is $t_x^2 \hat V(\hat B)$
\item We often use ratio estimators for means, letting the
auxiliary variable $x$ be an indicator (1 or 0) variable for whether
or not unit $i$ is in the sample\\
----Here, $\hat t_x$ is an estimator of the population size and the
ratio is the estimator of the population total divided by the
estimated population size
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Stratification and Clustering}
\bi
\item $H$ strata, within stratum $h$, select $n_h$ psu's
\item stratification: $\hat t = \sum\limits_{h=1}^H \hat t_h$
\item clustering part: $\hat t_h = \sum\limits_{i\in S_h}
\ds\frac{\hat t_{hi}}{\pi_{hi}}$
\item final total $\hat t = \sum\limits_{h=1}^H \sum\limits_{j\in S_h} \ds\frac{\hat
t_{hi}}{\pi_{hi}}$
\item $\hat V(\hat t) = \sum\limits_{h=1}^H \hat V(\hat t_h)$
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Any design}

Let $w_i$ be weight for unit $i$
\bi
\item $\hat t = \sum\limits_{i\in S} w_i y_i$
\item $\hat {\bar y} = \sum\limits_{i\in S} w_i
y_i/\sum\limits_{i\in S} w_i$\\
-----$\sum\limits_{i\in S}w_i$ is used to estimate the population size,
generally, not equal to true population size, except SRS, but should
be close.
\ei
\end{frame}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}

{\bf Example 7.1:} Malaria in The Gambia\\
Malaria is a serious health problem in The Gambia. Malaria morbidity
can be reduced by using bed nets that are impregnated with
insecticide, but this is only effective if the bed nets are in
widespread use. In 1991, a nationwide survey was designed to
estimate the prevalence of bed net use in rural areas. The survey is
described and results are reported in D'Alessandro et al.(1994).\\
\begin{tabular}{ccc}
Stage & Stratification & Sampling units \\
\hline
1 & Region & District  \\
2 & PHC/non-PHC & Village \\
3 & & Compound\\
\hline
\end{tabular}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\bi
\item The sampling frame consists: all rural villages of fewer than
3000 people in The Gambia.
\item The villages were stratified by three geographic regions
(eastern, central, and western).
\item In each region (eastern,
central, western), four districts were chosen with probability
proportional to the district population (used the 1983 census).
\item The district were stratified by whether the village had a
public health clinic (PHC) or not. Four villages were chosen, with
probability proportional to the 1983 census population: two PHC
villages and two non-PHC villages.
\item Six compounds were chosen
randomly from each village, and a researcher recorded the number of
beds and nets, along with other information, for each compound.
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}

{\bf Simplicity in Survey Design}
\bi
\item
 A simpler design giving the same amount
of information per dollar spent is almost always to be preferred to
a more complicated design\\
----often easier to administer and analyze, and data from the
simple designed survey are less likely to be analyzed incorrectly by
subsequent analysts.
\item Choose a more complex design if it is really more efficient and practical.
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Sampling Weights}
\bi
\item Stratified random sampling
$$\hat t_{\mbox{str}} = \sum\limits_{h=1}^H \sum\limits_{j \in S_h}
w_{hj} y_{hj}$$ $\pi_{hj} = n_h /N_h$: probability of selecting the
$j$th unit in the $h$th stratum to be in the sample\\
$w_{hj} = N_h/n_h$
$$\bar y_{\mbox{str}} = \ds\frac {\sum\limits_{h=1}^H \sum\limits_{j \in S_h}
w_{hj} y_{hj}} {\sum\limits_{h=1}^H \sum\limits_{j \in S_h}
w_{hj}}$$
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\bi
\item Cluster sampling with equal probabilities {\small
$$w_{ij} =\ds\frac{N M_i}{n m_i}= \ds\frac{1}{\mbox{probability that the jth ssu in the ith psu is in
the sample}}.$$
$$\hat t = \sum\limits_{i \in S} \sum\limits_{j\in S_i} w_{ij}
y_{ij}$$
$$\bar y = \ds\frac {\hat t}{\sum\limits_{i \in S} \sum\limits_{j\in S_i}
w_{ij}}$$

\item Cluster sampling with unequal probabilities}\\
Let $w_i$ be the weight that the $i$th psu is in the sample\\
and $w_{j|i}$ be the weight that the $j$th ssu is in the
sample given that the $i$th psu is in the sample\\
 the sampling weights are $$w_{ij} =w_i\times w_{j|i}$$
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\bi
\item Three-stage cluster sampling\\
$w_p$: weight for the psu\\
$w_{s|p}$: weight for the ssu\\
$w_{t|s,p}$: weight associated with the tsu (tertiary sampling
unit)\\
the overall sampling weight for an observation unit is $$w = w_p
\times w_{s|p} \times w_{t|s,p}$$

\item Any design\\
$w_i$: weight for unit $i$
$$\hat t = \sum\limits_{i\in S} w_i y_i$$
 $$\hat {\bar y} = \sum\limits_{i\in S} w_i
y_i/\sum\limits_{i\in S} w_i$$
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
Comments:
\bi
\item All the information needed to construct point
estimators is contained in the sampling weights
\item However, the
sampling weights give no information on how to find standard errors
of the estimators, and thus knowing the sampling weights alone will
not allow you to do inferential statistics
\item  Variances of estimators depend on the probabilities that any
pair of observation units is selected to be in the sample, and
requires more knowledge of the sampling design than given by weights
alone
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
{\bf Estimating a Distribution Function}

Review: For a discrete random variable $X$ with possible values
$x_1, x_2, \cdots, x_n$, a probability mass function is a function
such that
$$f(x_i) \ge 0, \quad \sum\limits_{i=1}^n f(x_i) = 1, \quad f(x_i) = P(X=x_i)$$


Cumulative distribution function is defined as
$$F(X) = P(X\le x) = \sum\limits_{x_i\le x} f(x_i)$$

Continuous case
$$F(x) = P(X\le x) = \int_{-\infty}^x f(x) dx$$
$$u =\int_{-\infty}^{+\infty} xdF(x) = \int_{-\infty}^{+\infty}
xf(x)dx$$
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
For sampling, suppose the values for the entire population of $N$
units are known
\bi
\item Probability mass function (pmf) of the finite population is
$$f(y) =\ds\frac{\mbox{number of units whose value is y}}{N}$$
\item Cumulative distribution function (cdf) of the finite
population is
$$F(y) =\ds\frac{\mbox{number of units with value} \le y}{N} =
\sum\limits_{x\le y} f(x)$$
\item Finite population mean $\bar y_U =\sum\limits_{\mbox{values of y in population}} yf(y)$
\ei
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\bi
\item Population median: any value $m$ with
$$F(m) \ge 1/2\quad \mbox{and} \quad P(Y\ge m) \ge 1/2$$
\item In general, $x$ is a 100rth percentile if
$$F(x) \ge r \quad\mbox{and} \quad P(Y\ge x) \ge 1-r$$
\item Population variance
\begin{eqnarray*}
S^2 &=& \ds\frac 1 {N-1} \sum\limits_{i=1}^N (y_i -\bar y_U)^2\\
&=& \ds\frac N {N-1} \sum\limits_y f(y) \left[y-\sum\limits_{y} y
f(y)\right]^2\\
&=& \ds\frac{N}{N-1} \left[\sum\limits_y y^2 f(y) -
\left(\sum\limits_y y f(y)\right)^2\right] \end{eqnarray*} 
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}

{\bf Estimators}
\bi
\item Empirical probability mass function (epmf)
$$\hat f(y) = \ds\frac{\sum\limits_{i\in S,y_i =y} w_i}{\sum\limits_{i\in
S}w_i}$$
\item Empirical distribution function (ecdf)
$$\hat F(y) = \sum\limits_{x \le y} \hat f(x)$$
or $$\hat F(y) = \ds\frac{\sum\limits_{i\in S, y_i \le y}
w_i}{\sum\limits_{i\in S}w_i}$$
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\bi
\item Estimate the median by $\hat F(y) \approx 1/2$
\item Estimate $\bar y_U$ by
$$\hat {\bar y} = \sum y\hat f(y) = \ds\frac{\sum\limits_{i\in S} y_i w_i}{\sum\limits_{i\in S} w_i}$$
\item Estimate population variance by
$$s^2 =\ds\frac{N}{N-1} \left[\sum\limits_y y^2 \hat f(y) - \left(\sum\limits_y y
\hat f(y)\right)^2\right]$$
\ei
\end{frame}

\begin{frame}
{\bf Plotting Data from Complex Survey}
\bi
\item plots commonly used for SRS can mislead when applied to raw
data from non-self-weighting samples
\item clustering causes numerous difficulties in plotting data from a
complex survey, because the clustering structure as well as possible
unequal weighting must be displayed in the graphs
\item the problems
are compounded because data sets from surveys are often very large
and involve several layers of clustering
\item plot data both with and without weights to see the effect of the
weights
\item plot data separately for each stratum and for each psu,
if possible, to examine variability in the responses
\item incorporating the weights into the graphics
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}

{\bf Example:}

Consider an artificial population of 1000 men and 1000 women stored in the file \texttt{htpop.dat}. Each individual's height is measured to the nearest centimeter (cm). 

Two samples are drawn from this population:
\bi
\item \texttt{htsrs.dat}: A simple random sample (SRS) of size 200.
\item \texttt{htstrat.dat}: A stratified random sample consisting of 160 women and 40 men.
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# Empirical pmf for stratified sample of heights}
\hlcom{# define sampling weight}
\hlkwd{library}\hlstd{(survey)}
\end{alltt}


{\ttfamily\noindent\itshape\color{messagecolor}{\#\# Loading required package: grid}}

{\ttfamily\noindent\itshape\color{messagecolor}{\#\# Loading required package: Matrix}}

{\ttfamily\noindent\itshape\color{messagecolor}{\#\# Loading required package: survival}}

{\ttfamily\noindent\itshape\color{messagecolor}{\#\# \\\#\# Attaching package: 'survey'}}

{\ttfamily\noindent\itshape\color{messagecolor}{\#\# The following object is masked from 'package:graphics':\\\#\# \\\#\# \ \ \ \ dotchart}}\begin{alltt}
\hlkwd{library}\hlstd{(SDAResources)}
\hlkwd{library}\hlstd{(tibble)}
\hlstd{htstrat}\hlopt{$}\hlstd{sampwt} \hlkwb{<-} \hlnum{1000}\hlopt{/}\hlkwd{sum}\hlstd{(htstrat}\hlopt{$}\hlstd{gender}\hlopt{==}\hlstr{"F"}\hlstd{)}
\hlstd{htstrat}\hlopt{$}\hlstd{sampwt[htstrat}\hlopt{$}\hlstd{gender}\hlopt{==}\hlstr{"M"}\hlstd{]} \hlkwb{<-}
  \hlnum{1000}\hlopt{/}\hlkwd{sum}\hlstd{(htstrat}\hlopt{$}\hlstd{gender}\hlopt{==}\hlstr{"M"}\hlstd{)}
\hlkwd{head}\hlstd{(htstrat)}
\end{alltt}
\begin{verbatim}
## # A tibble: 6 x 4
##      rn height gender sampwt
##   <dbl>  <dbl> <chr>   <dbl>
## 1   201    166 F        6.25
## 2   965    163 F        6.25
## 3   490    166 F        6.25
## 4   249    155 F        6.25
## 5   260    154 F        6.25
## 6   324    160 F        6.25
\end{verbatim}
\begin{alltt}
\hlcom{# use function emppmf to calculate pmf}
\hlstd{strresult} \hlkwb{<-} \hlkwd{emppmf}\hlstd{(htstrat}\hlopt{$}\hlstd{height,htstrat}\hlopt{$}\hlstd{sampwt)}
\hlkwd{print}\hlstd{(}\hlkwd{as_tibble}\hlstd{(}\hlkwd{head}\hlstd{(strresult)),} \hlkwc{n} \hlstd{=} \hlnum{6}\hlstd{)}
\end{alltt}
\begin{verbatim}
## # A tibble: 45 x 2
##    vals      epmf
##   <dbl> <dbl[1d]>
## 1   142   0.00312
## 2   144   0.00312
## 3   147   0.00312
## 4   148   0.00625
## 5   149   0.00312
## 6   150   0.00625
## # i 39 more rows
\end{verbatim}
\begin{alltt}
\hlcom{# plot }
\hlkwd{par}\hlstd{(}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{)}
\hlkwd{plot}\hlstd{(strresult}\hlopt{$}\hlstd{vals, strresult}\hlopt{$}\hlstd{epmf,}\hlkwc{type}\hlstd{=}\hlstr{"h"}\hlstd{,}
 \hlkwc{xlab}\hlstd{=}\hlstr{"Height Value, y (cm)"}\hlstd{,}
 \hlkwc{ylab}\hlstd{=}\hlstr{"Empirical pdf"}\hlstd{,}\hlkwc{lwd}\hlstd{=}\hlnum{1.2}\hlstd{,}
 \hlkwc{main}\hlstd{=}\hlstr{"Empirical pdf for stratified sample of heights (weighted)"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/pdf-1} 
\end{knitrout}
\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# Empirical cdf of height for data htpop, and for data htstrat, }
\hlcom{# with and without weights}
\hlstd{d0710} \hlkwb{<-} \hlkwd{svydesign}\hlstd{(}\hlkwc{id} \hlstd{=} \hlopt{~}\hlnum{1}\hlstd{,} \hlkwc{strata} \hlstd{=} \hlopt{~}\hlstd{gender,}
                   \hlkwc{fpc} \hlstd{=} \hlkwd{c}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{1000}\hlstd{,}\hlnum{160}\hlstd{),}\hlkwd{rep}\hlstd{(}\hlnum{1000}\hlstd{,}\hlnum{40}\hlstd{)),}
                   \hlkwc{data} \hlstd{= htstrat)}
\hlstd{cdf.weighted}\hlkwb{<-}\hlkwd{svycdf}\hlstd{(}\hlopt{~}\hlstd{height, d0710)}
\hlstd{cdf.weighted}
\end{alltt}
\begin{verbatim}
## Weighted ECDFs: svycdf(~height, d0710)
\end{verbatim}
\begin{alltt}
\hlcom{## evaluate the function for height 144}
\hlstd{cdf.weighted[[}\hlnum{1}\hlstd{]](}\hlnum{144}\hlstd{)}
\end{alltt}
\begin{verbatim}
## [1] 0.00625
\end{verbatim}
\begin{alltt}
\hlcom{## compare to population and unweighted sample ecdfs.}
\hlstd{cdf.pop}\hlkwb{<-}\hlkwd{ecdf}\hlstd{(htpop}\hlopt{$}\hlstd{height)}    \hlcom{# ecdf for population}
\hlstd{cdf.samp}\hlkwb{<-}\hlkwd{ecdf}\hlstd{(htstrat}\hlopt{$}\hlstd{height)} \hlcom{# unweighted ecdf of sample}
\hlkwd{par}\hlstd{(}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{,}\hlkwc{mar}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{5.1}\hlstd{,}\hlnum{4.1}\hlstd{,}\hlnum{2.1}\hlstd{,}\hlnum{2.1}\hlstd{))}
\hlkwd{plot}\hlstd{(cdf.pop,} \hlkwc{do.points} \hlstd{=} \hlnum{FALSE}\hlstd{,}
     \hlkwc{xlab}\hlstd{=}\hlstr{"Height value y"}\hlstd{,}\hlkwc{ylab}\hlstd{=}\hlstr{"Empirical cdf"}\hlstd{,}\hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{135}\hlstd{,}\hlnum{205}\hlstd{),}
     \hlkwc{lwd}\hlstd{=}\hlnum{2.5}\hlstd{,}
     \hlkwc{main}\hlstd{=}\hlstr{"Empirical cdfs for population and sample"}\hlstd{)}
\hlkwd{lines}\hlstd{(cdf.samp,} \hlkwc{col}\hlstd{=}\hlstr{"red"}\hlstd{,} \hlkwc{do.points} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{lwd}\hlstd{=}\hlnum{1.8}\hlstd{)}
\hlkwd{lines}\hlstd{(cdf.weighted[[}\hlnum{1}\hlstd{]],} \hlkwc{do.points} \hlstd{=} \hlnum{FALSE}\hlstd{,}  \hlkwc{col} \hlstd{=}\hlstr{"purple"}\hlstd{,}\hlkwc{lwd}\hlstd{=}\hlnum{1}\hlstd{)}
\hlkwd{legend}\hlstd{(}\hlstr{"topleft"}\hlstd{,} \hlkwc{legend}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"Population"}\hlstd{,} \hlstr{"Sample unweighted"}\hlstd{,}
    \hlstr{"Sample weighted"}\hlstd{),}\hlkwc{col}\hlstd{=}\hlkwd{c}\hlstd{(}\hlstr{"black"}\hlstd{,} \hlstr{"red"}\hlstd{,} \hlstr{"purple"}\hlstd{),}
    \hlkwc{lwd}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{2.5}\hlstd{,}\hlnum{1.8}\hlstd{,}\hlnum{1}\hlstd{),}\hlkwc{cex}\hlstd{=}\hlnum{0.8}\hlstd{,}\hlkwc{bty}\hlstd{=}\hlstr{"n"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/cdf-1} 
\end{knitrout}
\begin{frame}

{\bf Examples 7.5, 7.6, and 9.12}
\begin{table}
\caption{Estimates from population and different samples}
\begin{tabular}{ccccc} \hline Quantity & Population & SRS &
Stratified & Stratified\\
& && no weights & with weights\\
\hline Mean & 168.6 & 168.9 & 164.6 & 169.0\\
Median & 168 & 169 & 163 & 168 \\
25th percentile & 160 & 160& 157 &161\\
90th percentile & 184 & 184 & 178 & 182\\
variance & 124.5 &122.6 &93.4 &116.8 \\
\hline
\end{tabular}
\end{table}
\end{frame}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{##### Example 7.5}

\hlkwd{data}\hlstd{(htsrs)}
\hlstd{dhtsrs}\hlkwb{<-}\hlkwd{svydesign}\hlstd{(}\hlkwc{id} \hlstd{=} \hlopt{~}\hlnum{1}\hlstd{,}\hlkwc{weights}\hlstd{=}\hlkwd{rep}\hlstd{(}\hlnum{2000}\hlopt{/}\hlnum{200}\hlstd{,}\hlnum{200}\hlstd{),}
    \hlkwc{fpc}\hlstd{=}\hlkwd{rep}\hlstd{(}\hlnum{2000}\hlstd{,}\hlnum{200}\hlstd{),} \hlkwc{data}\hlstd{=htsrs)}
\hlcom{# cdf treated as step function, gives values in Table 7.1 of SDA}
\hlkwd{svyquantile}\hlstd{(}\hlopt{~}\hlstd{height, dhtsrs,} \hlkwc{quantiles}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0.25}\hlstd{,}\hlnum{0.5}\hlstd{,}\hlnum{0.75}\hlstd{,}\hlnum{0.9}\hlstd{),}
            \hlkwc{ties} \hlstd{=} \hlstr{"discrete"}\hlstd{)}
\end{alltt}
\begin{verbatim}
## $height
##      quantile ci.2.5 ci.97.5        se
## 0.25      160    159     163 1.0142211
## 0.5       169    168     171 0.7606658
## 0.75      176    174     179 1.2677764
## 0.9       184    182     187 1.2677764
## 
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"
\end{verbatim}
\begin{alltt}
\hlcom{# interpolated quantiles (usually preferred method)}
\hlkwd{svyquantile}\hlstd{(}\hlopt{~}\hlstd{height, dhtsrs,} \hlkwc{quantiles}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0.25}\hlstd{,}\hlnum{0.5}\hlstd{,}\hlnum{0.75}\hlstd{,}\hlnum{0.9}\hlstd{),}
            \hlkwc{ties} \hlstd{=} \hlstr{"rounded"}\hlstd{)}
\end{alltt}
\begin{verbatim}
## $height
##      quantile ci.2.5 ci.97.5        se
## 0.25      160    159     163 1.0142211
## 0.5       169    168     171 0.7606658
## 0.75      176    174     179 1.2677764
## 0.9       184    182     187 1.2677764
## 
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"
\end{verbatim}
\end{kframe}
\end{knitrout}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{##### Examples 7.6 and 9.12}

\hlkwd{data}\hlstd{(htstrat)}
\hlstd{popsize_recode} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{'F'} \hlstd{=} \hlnum{1000}\hlstd{,} \hlstr{'M'} \hlstd{=} \hlnum{1000}\hlstd{)}
\hlcom{# create a new variable popsize for population size}
\hlstd{htstrat}\hlopt{$}\hlstd{popsize}\hlkwb{<-}\hlstd{popsize_recode[htstrat}\hlopt{$}\hlstd{gender]}
\hlkwd{head}\hlstd{(}\hlkwd{as.data.frame}\hlstd{(htstrat))}
\end{alltt}
\begin{verbatim}
##    rn height gender popsize
## 1 201    166      F    1000
## 2 965    163      F    1000
## 3 490    166      F    1000
## 4 249    155      F    1000
## 5 260    154      F    1000
## 6 324    160      F    1000
\end{verbatim}
\begin{alltt}
\hlcom{# design object}
\hlcom{# svydesign calculates the weights here from the fpc argument}
\hlstd{dhtstrat}\hlkwb{<-}\hlkwd{svydesign}\hlstd{(}\hlkwc{id} \hlstd{=} \hlopt{~}\hlnum{1}\hlstd{,} \hlkwc{strata} \hlstd{=} \hlopt{~}\hlstd{gender,}
    \hlkwc{fpc} \hlstd{=} \hlopt{~}\hlstd{popsize,}\hlkwc{data} \hlstd{= htstrat)}
\hlcom{# ties = "discrete" gives values in Table 7.1 of SDA}
\hlkwd{svyquantile}\hlstd{(}\hlopt{~}\hlstd{height, dhtstrat,} \hlkwd{c}\hlstd{(}\hlnum{0.25}\hlstd{,}\hlnum{0.5}\hlstd{,}\hlnum{0.75}\hlstd{,}\hlnum{0.9}\hlstd{),}
            \hlkwc{ties} \hlstd{=} \hlstr{"discrete"}\hlstd{)}
\end{alltt}
\begin{verbatim}
## $height
##      quantile ci.2.5 ci.97.5        se
## 0.25      161    160     163 0.7606423
## 0.5       168    166     171 1.2677372
## 0.75      177    174     179 1.2677372
## 0.9       182    180     191 2.7890219
## 
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"
\end{verbatim}
\begin{alltt}
\hlcom{# ties = "rounded" gives values in Example 9.12 of SDA}
\hlkwd{svyquantile}\hlstd{(}\hlopt{~}\hlstd{height, dhtstrat,} \hlkwd{c}\hlstd{(}\hlnum{0.25}\hlstd{,}\hlnum{0.5}\hlstd{,}\hlnum{0.75}\hlstd{,}\hlnum{0.9}\hlstd{),}
            \hlkwc{ties} \hlstd{=} \hlstr{"rounded"}\hlstd{,}
            \hlkwc{ci}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{interval.type} \hlstd{=} \hlstr{"beta"}\hlstd{)}
\end{alltt}
\begin{verbatim}
## $height
##      quantile ci.2.5 ci.97.5        se
## 0.25      161    160     163 0.7606423
## 0.5       168    166     171 1.2677372
## 0.75      177    174     179 1.2677372
## 0.9       182    179     190 2.7890219
## 
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"
\end{verbatim}
\end{kframe}
\end{knitrout}

\begin{frame}

{\bf Examples 7.6 and 9.12}
\bi
\item data: htstrat
\item Histogram with smoothed density function
\item Boxplots
\ei
\end{frame}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{##### Examples 7.6 and 9.12}
\hlcom{# histogram and smoothed density function}
\hlkwd{data}\hlstd{(htstrat)}
\hlcom{# set graphics parameters, 1*2 plots, axis labels horizontal}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{),}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{,}\hlkwc{mar}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{2.1}\hlstd{,}\hlnum{4.1}\hlstd{,}\hlnum{2.1}\hlstd{,}\hlnum{0.3}\hlstd{))}
\hlcom{# Histogram overlaid with kernel density curve }
\hlcom{# (without weight information)}
\hlcom{# Displays the sample values, but does not estimate }
\hlcom{# population histogram}
\hlcom{# freq=FALSE changes the vertical axis to density}
\hlcom{# breaks tell how many breakpoints to use}
\hlkwd{hist}\hlstd{(htstrat}\hlopt{$}\hlstd{height,}\hlkwc{main}\hlstd{=}\hlstr{"Without weights"}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{"Height (cm)"}\hlstd{,}
     \hlkwc{breaks} \hlstd{=} \hlnum{10}\hlstd{,} \hlkwc{col}\hlstd{=}\hlstr{"gray90"}\hlstd{,} \hlkwc{freq}\hlstd{=}\hlnum{FALSE}\hlstd{,}
     \hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{140}\hlstd{,}\hlnum{200}\hlstd{),} \hlkwc{ylim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{0.045}\hlstd{))}
\hlcom{# overlaid with kernel density curve}
\hlkwd{lines}\hlstd{(}\hlkwd{density}\hlstd{(htstrat}\hlopt{$}\hlstd{height),}\hlkwc{lty}\hlstd{=}\hlnum{1}\hlstd{,}\hlkwc{lwd}\hlstd{=}\hlnum{2}\hlstd{)}
\hlcom{# Histogram (with weight information)}
\hlcom{# create survey design object, weights calculated from fpc here}
\hlstd{d0710} \hlkwb{<-} \hlkwd{svydesign}\hlstd{(}\hlkwc{id} \hlstd{=} \hlopt{~}\hlnum{1}\hlstd{,} \hlkwc{strata} \hlstd{=} \hlopt{~}\hlstd{gender,}
                   \hlkwc{fpc} \hlstd{=} \hlkwd{c}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{1000}\hlstd{,}\hlnum{160}\hlstd{),}\hlkwd{rep}\hlstd{(}\hlnum{1000}\hlstd{,}\hlnum{40}\hlstd{)),}
                \hlkwc{data} \hlstd{= htstrat)}
\hlstd{d0710}
\end{alltt}
\begin{verbatim}
## Stratified Independent Sampling design
## svydesign(id = ~1, strata = ~gender, fpc = c(rep(1000, 160), 
##     rep(1000, 40)), data = htstrat)
\end{verbatim}
\begin{alltt}
\hlkwd{svyhist}\hlstd{(}\hlopt{~}\hlstd{height,d0710,} \hlkwc{main}\hlstd{=}\hlstr{"With weights"}\hlstd{,}\hlkwc{xlab} \hlstd{=} \hlstr{"Height (cm)"}\hlstd{,}
        \hlkwc{breaks} \hlstd{=} \hlnum{10}\hlstd{,} \hlkwc{col}\hlstd{=}\hlstr{"gray90"}\hlstd{,} \hlkwc{freq}\hlstd{=}\hlnum{FALSE}\hlstd{,}
        \hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{140}\hlstd{,}\hlnum{200}\hlstd{),} \hlkwc{ylim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{0.045}\hlstd{))}
\hlstd{dens1}\hlkwb{<-}\hlkwd{svysmooth}\hlstd{(}\hlopt{~}\hlstd{height,d0710,}\hlkwc{bandwidth}\hlstd{=}\hlnum{5}\hlstd{)}
\hlkwd{lines}\hlstd{(dens1,}\hlkwc{lwd}\hlstd{=}\hlnum{2}\hlstd{)} \hlcom{# draw the density line}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/hist_htstrat-1} 
\end{knitrout}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# boxplot}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{),}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{,}\hlkwc{mar}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{2.1}\hlstd{,}\hlnum{4.1}\hlstd{,}\hlnum{2.1}\hlstd{,}\hlnum{0.3}\hlstd{))}
\hlcom{# boxplot (with weight information)}
\hlkwd{svyboxplot}\hlstd{(height}\hlopt{~}\hlnum{1}\hlstd{,d0710,}\hlkwc{ylab}\hlstd{=}\hlstr{"Height"}\hlstd{,}\hlkwc{xlab}\hlstd{=}\hlstr{" "}\hlstd{,}
           \hlkwc{main}\hlstd{=}\hlstr{"Full sample"}\hlstd{)}
\hlkwd{svyboxplot}\hlstd{(height}\hlopt{~}\hlstd{gender,d0710,}\hlkwc{ylab}\hlstd{=}\hlstr{"Height"}\hlstd{,}\hlkwc{xlab}\hlstd{=}\hlstr{"Gender"}\hlstd{,}
           \hlkwc{main}\hlstd{=}\hlstr{"Separately by gender"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/boxplot_htstrat-1} 
\end{knitrout}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}

{\bf Example 7.9:} National Health and Nutrition Examination Survey (NHANES),
(Centers for Disease Control and Prevention, 2017). 
In this example, we look at statistics about body mass index (BMI,
variable bmxbmi ) for adults age 20 and over and plots. 
\bi
\item The \textit{weights} argument in \textit{svydesign} specifies the variable containing the final weights.\\\
 \textit{wtint2yr}: weight for the set of persons with interview data,\\
 \textit{wtmec2yr}: weight for the subset of interviewed persons who had a medical examination. \\
BMI is measured in the medical examination, so the appropriate weight variable to use is \textit{wtmec2yr}. 
\item The \textit{strata} and \textit{id} arguments are used as before, except now we include both of them. The \texttt{strata=\~{}sdmvstra} argument says that \textit{sdmvstra} is the variable giving the stratum membership. The psus specified in \texttt{id=\~{}sdmvpsu} are the first-stage sampling units.
\ei
\end{frame}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\bi
\item 
In the NHANES data, the two psus in each stratum are labeled as `1' and `2'. The \texttt{nest=TRUE} argument says that psu labels are nested within strata---that is, multiple strata have the same psu labels.  Typing \texttt{nest=TRUE} ensures that psu 1 in stratum 1 is recognized as being a different psu than psu 1 in stratum 2.

When there is one stratum, the results are the same if you have \texttt{nest=TRUE}, \texttt{nest=FALSE}, or simply omit the \textit{nest} argument.
\item
No \textit{fpc} argument is included in \textit{svydesign}. With complex samples such as NHANES, we usually want to calculate the with-replacement variance, which requires only psu-level information.
\ei
 
\end{frame}
\begin{knitrout}\scriptsize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{### Computing Estimates from Stratified Multistage Samples ####}
\hlcom{# note: load package in this order to get correct data nhanes}
\hlcom{# library(survey)}
\hlcom{# library(SDAResources)}
\hlkwd{data}\hlstd{(nhanes)}
\hlkwd{nrow}\hlstd{(nhanes)} \hlcom{#9971}
\end{alltt}
\begin{verbatim}
## [1] 9971
\end{verbatim}
\begin{alltt}
\hlkwd{names}\hlstd{(nhanes)}
\end{alltt}
\begin{verbatim}
##  [1] "sdmvstra" "sdmvpsu"  "wtint2yr" "wtmec2yr" "ridstatr" "ridageyr"
##  [7] "ridagemn" "riagendr" "ridreth3" "dmdeduc2" "dmdfmsiz" "indfmpir"
## [13] "bmxwt"    "bmxht"    "bmxbmi"   "bmxwaist" "bmxleg"   "bmxarml" 
## [19] "bmxarmc"  "bmdavsad" "lbxtc"    "bpxpls"   "sbp"      "dbp"     
## [25] "bpread"
\end{verbatim}
\begin{alltt}
\hlcom{# count number of observations with missing value}
\hlcom{# for ridageyr, bmxbmi}
\hlkwd{sum}\hlstd{(}\hlkwd{is.na}\hlstd{(nhanes}\hlopt{$}\hlstd{ridageyr))} \hlcom{# ridageyr gives age in years}
\end{alltt}
\begin{verbatim}
## [1] 0
\end{verbatim}
\begin{alltt}
\hlkwd{sum}\hlstd{(}\hlkwd{is.na}\hlstd{(nhanes}\hlopt{$}\hlstd{bmxbmi))}   \hlcom{# bmxbmi gives BMI}
\end{alltt}
\begin{verbatim}
## [1] 1215
\end{verbatim}
\end{kframe}
\end{knitrout}
\begin{frame}
\bi
\item
 \textit{subset} function: define domains. This specifies that estimates are desired for the domain of persons age 20 and older having data for BMI (with \textit{age20d}=1), and carries the stratification and clustering information from the full design over for analyzing the subset. 
 \item If you just created a subset of the data consisting of the observations having \textit{ridageyr} $\geq 20$, in some instances (for example, when some psus have no members of the domain), the standard errors would be incorrect; by using the \textit{subset} function,  the correct standard errors are calculated.
\item  We exclude adults with missing values of \textit{bmxbmi} from the domain of interest with \textit{age20d}=1. The estimates are computed from the adults who have data. If the domain contained missing values, we would need to include option \texttt{na.rm=TRUE} in the \textit{svymean} function to be able to calculate statistics.
\ei
\end{frame}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# define age20d}
\hlstd{nhanes}\hlopt{$}\hlstd{age20d}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwd{nrow}\hlstd{(nhanes))}
\hlstd{nhanes}\hlopt{$}\hlstd{age20d[nhanes}\hlopt{$}\hlstd{ridageyr} \hlopt{>=}\hlnum{20} \hlopt{& !}\hlkwd{is.na}\hlstd{(nhanes}\hlopt{$}\hlstd{bmxbmi)]}\hlkwb{<-}\hlnum{1}
\hlcom{# check missing value counts for new variables}
\hlkwd{sum}\hlstd{(}\hlkwd{is.na}\hlstd{(nhanes}\hlopt{$}\hlstd{age20d))}
\end{alltt}
\begin{verbatim}
## [1] 0
\end{verbatim}
\begin{alltt}
\hlkwd{sum}\hlstd{(nhanes}\hlopt{$}\hlstd{age20d)} \hlcom{# how many records in domain?}
\end{alltt}
\begin{verbatim}
## [1] 5406
\end{verbatim}
\begin{alltt}
\hlkwd{head}\hlstd{(nhanes[,} \hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{,} \hlnum{15}\hlstd{,} \hlnum{26}\hlstd{)])}
\end{alltt}
\begin{verbatim}
## # A tibble: 6 x 6
##   sdmvstra sdmvpsu wtint2yr wtmec2yr bmxbmi age20d
##      <dbl>   <dbl>    <dbl>    <dbl>  <dbl>  <dbl>
## 1      125       1  134671.  135630.   27.8      1
## 2      125       1   24329.   25282.   30.8      1
## 3      131       1   12400.   12576.   28.8      1
## 4      131       1  102718.  102079.   42.4      1
## 5      126       2   17628.   18235.   20.3      1
## 6      128       1   11252.   10879.   28.6      1
\end{verbatim}
\end{kframe}
\end{knitrout}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# stratified cluster design }
\hlstd{d0709} \hlkwb{<-} \hlkwd{svydesign}\hlstd{(}\hlkwc{id} \hlstd{=} \hlopt{~}\hlstd{sdmvpsu,} \hlkwc{strata}\hlstd{=}\hlopt{~}\hlstd{sdmvstra,}
  \hlkwc{weights}\hlstd{=}\hlopt{~}\hlstd{wtmec2yr,} \hlkwc{nest}\hlstd{=}\hlnum{TRUE}\hlstd{,} \hlkwc{data} \hlstd{= nhanes)}
\hlcom{# domain estimation, age20+ }
\hlstd{d0709sub}\hlkwb{<-}\hlkwd{subset}\hlstd{(d0709, age20d} \hlopt{==}\hlnum{1}\hlstd{)}
\hlstd{d0709sub}
\end{alltt}
\begin{verbatim}
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (30) clusters.
## subset(d0709, age20d == 1)
\end{verbatim}
\begin{alltt}
\hlcom{# Request means and design effects}
\hlstd{nhmeans}\hlkwb{<-}\hlkwd{svymean}\hlstd{(}\hlopt{~}\hlstd{bmxbmi, d0709sub,} \hlkwc{deff}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\hlkwd{degf}\hlstd{(d0709sub)}
\end{alltt}
\begin{verbatim}
## [1] 15
\end{verbatim}
\begin{alltt}
\hlstd{nhmeans}
\end{alltt}
\begin{verbatim}
##           mean      SE   DEff
## bmxbmi 29.3891  0.2532 7.1248
\end{verbatim}
\begin{alltt}
\hlkwd{confint}\hlstd{(nhmeans,}\hlkwc{df}\hlstd{=}\hlkwd{degf}\hlstd{(d0709sub))}
\end{alltt}
\begin{verbatim}
##           2.5 %   97.5 %
## bmxbmi 28.84942 29.92878
\end{verbatim}
\end{kframe}
\end{knitrout}

\begin{frame}

The mean BMI for adults aged 20 and over is **29.389**, calculated using the design object \textit{d0709sub}, with a 95\% confidence interval of [28.849, 29.929].

\vspace{0.3in}

Next, we present bubble plots along with a smoothed trend line. The area of each circle in the plot is proportional to the sum of the weights for the sample observations.
\end{frame}


\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{### Bubble plots}
\hlkwd{par}\hlstd{(}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{)} \hlcom{# make tick mark labels horizontal}
\hlkwd{svyplot}\hlstd{(bmxbmi}\hlopt{~}\hlstd{ridageyr,} \hlkwc{design}\hlstd{=d0709,} \hlkwc{style}\hlstd{=}\hlstr{"bubble"}\hlstd{,} \hlkwc{inches}\hlstd{=}\hlnum{0.03}\hlstd{,}
        \hlkwc{xlab}\hlstd{=}\hlstr{"Age(years)"}\hlstd{,}\hlkwc{ylab}\hlstd{=}\hlstr{"Body Mass Index"}\hlstd{,}
        \hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{80}\hlstd{),}\hlkwc{ylim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{10}\hlstd{,}\hlnum{70}\hlstd{),}
        \hlkwc{main}\hlstd{=}\hlstr{"Weighted bubble plot of BMI versus age"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/nhanes3-1} 
\end{knitrout}

\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{### Plot data for a domain}
\hlcom{# define subset}
\hlstd{d0709subA}\hlkwb{<-}\hlkwd{subset}\hlstd{(d0709, ridreth3}\hlopt{==}\hlnum{6}\hlstd{)}
\hlkwd{par}\hlstd{(}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{)} \hlcom{# make tick mark labels horizontal}
\hlkwd{svyplot}\hlstd{(bmxbmi}\hlopt{~}\hlstd{ridageyr,} \hlkwc{design}\hlstd{=d0709subA,}
        \hlkwc{style}\hlstd{=}\hlstr{"bubble"}\hlstd{,}\hlkwc{inches} \hlstd{=} \hlnum{0.03}\hlstd{,}
        \hlkwc{xlab}\hlstd{=}\hlstr{"Age(years)"}\hlstd{,}\hlkwc{ylab}\hlstd{=}\hlstr{"Body Mass Index"}\hlstd{,}
        \hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{80}\hlstd{),}\hlkwc{ylim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{10}\hlstd{,}\hlnum{70}\hlstd{),}
\hlkwc{main}\hlstd{=}\hlstr{"Weighted bubble plot of BMI versus age for Asian Americans"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/nhanes4-1} 
\end{knitrout}
\begin{knitrout}\footnotesize
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{### Smoothed trend line for mean}
\hlcom{# Smoothed trend line with bubble plot of BMI versus age}
\hlcom{# plot data bmxbmi~ridageyr}
\hlkwd{par}\hlstd{(}\hlkwc{las}\hlstd{=}\hlnum{1}\hlstd{)} \hlcom{# make tick mark labels horizontal}
\hlkwd{svyplot}\hlstd{(bmxbmi}\hlopt{~}\hlstd{ridageyr,} \hlkwc{design}\hlstd{=d0709,} \hlkwc{style}\hlstd{=}\hlstr{"bubble"}\hlstd{,}
        \hlkwc{basecol}\hlstd{=}\hlstr{"gray"}\hlstd{,}\hlkwc{inches}\hlstd{=}\hlnum{0.03}\hlstd{,}
        \hlkwc{xlab}\hlstd{=}\hlstr{"Age(years)"}\hlstd{,}\hlkwc{ylab}\hlstd{=}\hlstr{"Body Mass Index"}\hlstd{,}
        \hlkwc{xlim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{80}\hlstd{),}\hlkwc{ylim}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{10}\hlstd{,}\hlnum{70}\hlstd{),}
\hlkwc{main}\hlstd{=}\hlstr{"Smoothed trend line with bubble plot of BMI versus age"}\hlstd{)}
\hlcom{# plot smoothing trend line}
\hlcom{# library(KernSmooth)  # install and load the package}
\hlstd{smth}\hlkwb{<-}\hlkwd{svysmooth}\hlstd{(bmxbmi}\hlopt{~}\hlstd{ridageyr,d0709)}
\hlkwd{lines}\hlstd{(smth,}\hlkwc{lwd}\hlstd{=}\hlnum{2}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/nhanes5-1} 
\end{knitrout}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Disclaimer}

Slides are intended for a course based on the book: \textit{Sampling: Design and Analysis, Third Edition} by Sharon L. Lohr

All examples, datasets, and pages directly scanned and included in this chapter are from the textbook.

References to papers or books cited in these slides can be found in the Bibliography section of the textbook.
\end{frame}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
