## Load packages and data

```{r}
library(mase)
library(survey)
library(tidyverse)
data(api)
```

## Exploratory data analysis

```{r}
glimpse(apisrs)


ggplot(data = apisrs, 
       mapping = aes(x = meals, 
                     y = api00, 
                     color = stype)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "lm", se = FALSE)
```

## Estimation

### Horvitz-Thompson Estimator

```{r}

horvitzThompson(y = apisrs$api00,
          pi = 1/apisrs$pw, var_est = TRUE)
```

### Generalized Regression Estimator

```{r}
# Prepare the sample data
xsample <- select(apisrs, meals, stype) %>%
  as.data.frame()

# Grab the same auxiliary variables from the population dataset
xpop <- select(apipop, names(xsample))

# Fit the estimator
greg(y = apisrs$api00, xsample = xsample, xpop = xpop,
          pi = 1/apisrs$pw, var_est = TRUE, datatype = "raw")

# Look at the coefficient estimates
greg(y = apisrs$api00, xsample = xsample, xpop = xpop,
          pi = 1/apisrs$pw, var_est = TRUE,
     datatype = "raw")$coefficients

```

### GREG with the Elastic-Net

```{r}
gregElasticNet(y = apisrs$api00, xsample = xsample, xpop = xpop,
          pi = 1/apisrs$pw, var_est = TRUE, datatype = "raw")


gregElasticNet(y = apisrs$api00, xsample = xsample, xpop = xpop,
          pi = 1/apisrs$pw, var_est = TRUE,
          datatype = "raw")$coefficients
```

### GREG with a Regression Tree

```{r}
est <- gregTree(y = apisrs$api00, xsample = xsample, xpop = xpop,
          pi = 1/apisrs$pw, var_est = TRUE)
est

# Explore the tree
est$tree
```

### GREG with Summarized Auxiliary Data

```{r}
library(recipes)
xpop <- select(apipop, meals, stype)

# Assume you don't have access to xpop, just xpop_agg
xpop_agg <- recipe(~ meals + stype, data = xpop) %>%
  step_dummy(stype) %>%
  prep(training = xpop) %>%
  bake(new_data = NULL) %>%
  summarise_all(mean)
```

```{r}
# Wrangle the sample data
xsample <- recipe(~ meals + stype, data = apisrs) %>%
  step_dummy(stype) %>%
  prep(training = apisrs) %>%
  bake(new_data = NULL) %>%
  as.data.frame()

# Explore sample means
xsample %>%
  summarise_all(mean)

# Compute GREG
greg(y = apisrs$api00, xsample = xsample, xpop = xpop_agg,
          pi = 1/apisrs$pw, var_est = TRUE, datatype = "means")

```

### GREG with Logistic Regression

```{r}
apisrs <- mutate(apisrs, both_num = if_else(both == "Yes",  1, 0)) 

xsample <- select(apisrs, meals, stype) %>%
  as.data.frame()

greg(y = apisrs$both_num, xsample = xsample, xpop = xpop,
          pi = 1/apisrs$pw, var_est = TRUE, datatype = "raw",
     model = "logistic")
```

