\documentclass[12pt]{article}
\usepackage{verbatim, amssymb, amsmath, float, fullpage}
\begin{document}
\title{36-303 Sampling, Surveys \& Society \\ Homework 04 Solutions}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PROB1  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\setcounter{section}{1}

\section{Question 1 (Sarah Igo written lecture)}

(answers vary; I will post something more about this later --BJ)

\section{Question 2 (taken from S. Lohr's book)}
In this problem we have to think of the claims as our PSU. Each claim contains 215 fields that might have errors and $n=85$ {\it claims} have been sampled randomly from the population of $N=828$ claims.\\

The 85 claims have a total of $\sum_{i=1}^{85}t_i$ errors, where $t_i$ is the number of errors in claim $i$. An estimate of the total number of errors in our universe of $N=828$ claims is
\begin{eqnarray*}
\hat t = \frac{N}{n}=\frac{828}{85}\times 37 = 360.42,
\end{eqnarray*}
and an estimate of the error rate is
\begin{eqnarray*}
\hat{ \bar y} = {\hat t \over N \times 215} = {360 \over 828 \times 215} = 0.002025
\end{eqnarray*}
An estimate of the standard error of $\hat{\bar y}$ can be obtained from the formula
\begin{eqnarray*}
\hat{SE}[\hat{\bar y}] &=& \frac{1}{215}\sqrt{\left( 1- \frac{n}{N}\right)\frac{s^2_t}{n}}
\end{eqnarray*}
where $s_t^2$ is the sample variance of the psu's,
\begin{eqnarray*}
s_t^2 &=& \frac{1}{n-1}\sum_{i=1}^{n}\left( t_i - \frac{\hat t}{N}\right)^2\\
&=& \frac{1}{84}\sum_{i=1}^{84}\left( t_i - \frac{37}{85}\right)^2\\
&=& \frac{1}{84}\left[ 57\left(0-\frac{37}{85}\right)^2 + 
22\left(1-\frac{37}{85}\right)^2 +
4\left(2-\frac{37}{85}\right)^2 +
1\left(3-\frac{37}{85}\right)^2 +
1\left(4-\frac{37}{85}\right)^2
\right]\\
&=& 0.558263
\end{eqnarray*} 
Then, the SE for the estimate of the error rate is
\begin{eqnarray*}
\hat{SE}[\hat{\bar y}] &=& \frac{1}{215}\sqrt{\left( 1- \frac{85}{828}\right)\frac{0.558263}{85}} = 0.000357
\end{eqnarray*} 
\subsection{b)}
>From the previous part we have that $$\hat t = 360.42$$\\
An estimate of the standard error of this estimate ca be obtained using the formula
\begin{eqnarray*}
\hat{SE}[\hat{\bar y}] &=& 828\sqrt{\left( 1- 
\frac{85}{828}\right)\frac{s^2_t}{85}} = 63.565
\end{eqnarray*}

\subsection{c)}
In this case we have to think our universe as composed from `fields'. We have $N= 828 \times 215$ fields and our sample is composed from $n=85 \times 215$ fields. The error rate is obviously the same as in the previous case,
\begin{eqnarray}
\hat p_{SRS} = \frac{\mbox{fields with errors}}{n} = \frac{37}{85 \times 215} = 0.002025
\end{eqnarray}
The interesting thing happens when we compute the variance assuming that the sample is effectively a SRS {\it of fields}:
\begin{eqnarray*}
\hat V[\hat p_{SRS}] = \left( 1 - {85 \times 215 \over 828 \times 215 }\right) \frac{\hat p_{SRS}(1 - \hat p_{SRS})}{85 \times 215} = 9.92 \times 10^{-8}
\end{eqnarray*}

If we compare this estimate with the estimate obtained in part a), assuming cluster sampling,
\begin{eqnarray}
\hat V[ \hat p_{clus} ] = \left( \hat {SE}[\hat {\bar y}] \right)^2 = 0.000357^2= 1.27449 \times 10 ^{-7}
\end{eqnarray}
we see that this last variance estimate is bigger than the one computed assuming SRS. This is a general phenomenon when we have clustered samples. To achieve the same error levels, a clustered sample must be bigger than a SRS. This example also illustrate the problems of analyzing clustered samples using SRS methods: the SE of the estimates will be underestimated. This is dangerous because we  (and others) will think that our point estimates are better than they really are. 

\section{Question 3}
Creating the dataset in R:
\begin{verbatim}
strata <- data.frame(expand.grid(Sex=factor(c('M','F')), 
     College=factor(c('Eng','Lib'))),
  n_h = c(8,4,2,6),
  N_h = c(617,450,380,551),
  sam_w = NA,
  Pop_W = NA
)
strata$Pop_W <- strata$N_h / sum(strata$N_h)
strata$sam_w  <- strata$n_h / sum(strata$n_h)
HrsWk <- c(28,29,23,35,29,30,34,31,30,31,36,33,27,28,29,30,28,28,32,30)
data <- cbind(strata[rep(1:NROW(strata), strata$n_h),],HrsWk)
\end{verbatim}

\subsection{a)}
The mean of `Hrs/Wk' is
\begin{verbatim}
> mean(data$HrsWk)
[1] 30.05
\end{verbatim}
Since this is SRS without replacement, an estimate of the standard error of the mean is 
\begin{eqnarray*}
\hat {SE}[\bar y] &=& \sqrt{\left( 1 - {n\over N}\right)\frac{s^2}{n}}
\end{eqnarray*}
where $s^2$ is the sample variance. Computing this,
\begin{verbatim}
> fpc <- 1 - sum(strata$n_h)/sum(strata$N_h)
> sqrt(fpc) * sqrt(var(data$HrsWk) / NROW(data))
[1] 0.6634416
\end{verbatim}

\subsection{b)}
Computing the post-stratification weights,
\begin{verbatim}
data <- cbind(data,PSW = data$Pop_W/data$sam_w)
\end{verbatim}
And the weighted mean using the post-stratification weights is
\begin{verbatim}
> weighted.mean(data$HrsWk, w = data$PSW)
[1] 29.9111
\end{verbatim}
which is slightly lower than the one without using the population-level information.

\subsection{c)}
To estimate the SE using a first order Taylor series approximation we use the R procedure given in class (don't forget the finite population correction)
\begin{verbatim}
>tsv <- ts.variance(data$HrsWk, w = data$PSW)
> se.ts <- sqrt(fpc)*sqrt(tsv$var.ts)
> se.ts
[1] 0.6709889
\end{verbatim}

\section{Question 4}
To estimate the SE using the Jackknife technique we use the R function given in class,
\begin{verbatim}
> stacked_strata <- data.frame(stratum = paste(data$Sex, 
    data$College, sep='.'),data$N_h)
> jkv <- jk.variance(data$HrsWk, stacked_strata$stratum,
    unique(stacked_strata$data.N_h))
> jkv
$ybar.weighted
[1] 29.9111

$ybar.reps
[1] 29.9111

$var.jk
[1] 0.3419697
\end{verbatim}
Then, the jacknife estimate of the SE is
\begin{verbatim}
> se.jk <- sqrt(fpc)*sqrt(jkv$var.jk)
> se.jk
[1] 0.5818476
\end{verbatim}
\end{document}
