---
title: "Marginal Mean Weighting throught Stratification to Estimate the ATE"
output:
  html_document:
    df_print: paged
editor_options:
  chunk_output_type: inline
---
# Preparation Steps

## Load library that generates different types of weights based on the propensity score
```{r}
library(WeightIt) 
```

## Load data

```{r}
load(file="Chapter4_ssocs08_with_propensity_scores.Rdata")
```

## Define propensiy score formula

- The propensity score formula contains covariates and dummy missing data indicators.

```{r}
psFormula <- paste(covariateNames, collapse="+")
psFormula <- formula(paste("treat~",psFormula, sep=""))
print(psFormula)
```



# Obtain Marginal Mean Weights for the ATE

- weight it can estimate a variety of different types of propensity score weights
- weightit is able to estimate propensity scores with logistic regression and a variety of machine learning models
- weight it can incorporate sampling weights


```{r}
weights = weightit(formula=psFormula, #weighit is able to estimate propensity scores if covariates are provided 
                   data=SSOCS.data, #dataset
                   method="glm", #Use logistic regression to estimate propensity scores
                   estimand="ATE", #it can be ATE or ATT
                   stabilize = FALSE, #stabilizing weights can help reduce extreme weights 
                   s.weights="FINALWGT", #sampling weights
                   subclass=5) #specify that 5 strata are going to be used for MMWS
summary(weights)
                   
```

# Evaluate Covariate Balance

```{r}
library(cobalt) #package for evaluating covariate balance
balance.info = bal.tab(weights)
print(balance.info)

```

## Check if covariates met the WWC covariate balance criteria

### Check covariates with difference > 0.05
```{r}
balance.table = balance.info$Balance
print(row.names(balance.table)[abs(balance.table$Diff.Adj)>0.05]) 

```

### Check covariates with difference > 0.25

```{r}
print(row.names(balance.table)[abs(balance.table$Diff.Adj)>0.25]) 
```

# Estimate the average treatment effect

## Set up the survey design
```{r}
library(survey)
SSOCS.data$MMWS <- weights$weights #add weights to dataset
surveyDesignATE <- svydesign(ids=~0,strata=~STRATA,weights=~MMWS, data=SSOCS.data)
```

## Fit a linear regression model to proportions
```{r}
modelATE <- svyglm(percHarsh~treat, surveyDesignATE, family=gaussian())
summary(modelATE)
```


## Fit a logistic model  to proportions with quasibinomial family to allow for overdispersion

```{r}
modelATE2 <- svyglm(percHarsh~treat, surveyDesignATE, family=quasibinomial())
summary(modelATE2)
```

# Fit a model including the covariates that exceeded the 0.05 threshold

* Generate the outcome model formula
```{r}
covariatesToAdjust = row.names(balance.table)[abs(balance.table$Diff.Adj)>0.05]
yFormula <- paste(covariatesToAdjust, collapse="+")
yFormula <- formula(paste("percHarsh~treat+treat*(",yFormula,")", sep=""))
print(yFormula)
```
## Fit a logistic model  with covariates and interactions
*Testing interactions evaluates the homogeneity of regression slopes assumption of ANCOVA*

```{r}
yFormula = c("percHarsh ~ treat + treat * (pScores + C0190 + C0194 + 
    C0198 + C0202 + C0528 + C0538 + C0560 + C0562 + PROBWK08 + 
    FR_LVEL + C0542 + C0544 + C0556)")
modelATE3 <- svyglm(yFormula, surveyDesignATE, family=quasibinomial())
summary(modelATE3)
```

# Fit final model with treatment and covariates
*None of the interactions were significant so they were excluded from final model*

```{r}
yFormula = c("percHarsh ~ treat + pScores + C0190 + C0194 + 
    C0198 + C0202 + C0528 + C0538 + C0560 + C0562 + PROBWK08 + 
    FR_LVEL + C0542 + C0544 + C0556")
modelATE4 <- svyglm(yFormula, surveyDesignATE, family=quasibinomial())
summary(modelATE4)
```

-------------------------------------------

# PRACTICE

## Use a different method to estimate propensity scores 
- This can be done by changing the  method="glm" argument in the weightit function
*type help(weightit) to see methods available*

### Use generalized boosted modeling (gbm) to estimate propensity scores
```{r}

```

### Use the covariate balancing propensity score algorithm
```{r}


```


### Use the SuperLearner
```{r}


```