9: Machine Learning and Causal Inference

class: center, middle, inverse, title-slide

.title[
# 9: Machine Learning and Causal Inference
]
.subtitle[
## Quantitative Causal Inference
]
.author[
### <large>J. Seawright</large>
]
.institute[
### <small>Northwestern Political Science</small>
]
.date[
### May 29, 2025
]

---

class: center, middle

pre[class] {
  max-height: 200px;
}
</style>

---
### Heterogeneity

One of the central themes of the course has been causal heterogeneity. Can newer tools associated with machine learning contribute here?

---
### A (Quick) Historical Journey

`$$Y_{i} \sim \tau D_{i} + \beta \mathbb{X}_{i}$$`
where `$\tau = Y_{t} - Y_{c}$`.

---

`$$Y_{i} \sim \tau D_{i} + \beta f(\mathbb{X})_{i}$$`
where `$\tau = Y_{t} - Y_{c}$`.

---

Define:

`$$e(x) = E(W_{i} | \mathbb{X}_{i} = x)$$`
`$$m(x) = e(Y_{i} | \mathbb{X}_{i} = x) = f(x) + \tau e(x)$$`

---

Then we can rewrite our model as:

`$$Y_{i} - m(x) = \tau(W_{i} - e(x)) + \epsilon_{i}$$`
---

How do we get estimates of `$e$` and `$m$` that allow us to proceed with nonlinear equations without knowing the functional form?

We can run regressions, or we can just estimate them using nonparametric statistics or machine learning.

---

`$$Y_{i} \sim \tau(\mathbb{X}_{i}) D_{i} + \beta f(\mathbb{X})_{i}$$`

---

If we do what we were doing before *within neighborhoods* where `$\mathbb{X}_{i}$` is similar, then this will work. We will get multiple separate estimates that we can then analyze and pool.

Can we use random forests for this?

---
### CART

Optimally predicting `$Y$` based on `$\mathbb{X}$`, without assumptions of
additivity, linearity, etc.

---
### CART

1.  Start with the set of all cases, i.e., the root node.

2.  Search all values of each variable in `$\mathbb{X}$` for the binary
    split that maximizes homogeneity of `$Y$` for cases on each side of
    the split.

3.  If homogeneity of `$Y$` is sufficient or if the remaining sets of
    cases as the new final nodes are too small, stop. Otherwise, repeat
    step 2 for each of the current last-generation nodes.

---
### An Example

---

``` r
library(rpart)
```

---

``` r
dem.cart <- rpart(vdem_libdem ~ vdem_gender + vdem_corr + wdi_gdpcapcon2010 + wdi_mobile + wdi_gerp + une_surlgpef + wdi_fertility, data=qog_std_ts_jan22, na.action=na.omit)

plot(dem.cart)
text(dem.cart, use.n=TRUE)
```

---
### Random Forests

1.  Bootstrap the underlying data.

2.  Run CART, selecting a random subset of datapoints and variables at each
    decision node.

3.  Repeat several times, and find a way to average the results
    together.

---

``` r
library(randomForest)
```

```
## randomForest 4.7-1.2
```

```
## Type rfNews() to see new features/changes/bug fixes.
```

```
## 
## Attaching package: 'randomForest'
```

```
## The following object is masked from 'package:dplyr':
## 
##     combine
```

```
## The following object is masked from 'package:ggplot2':
## 
##     margin
```

---

``` r
dem.rf <- randomForest(vdem_libdem ~ wdi_gendeqr + bci_bci + wdi_gdpcapcon2010 + wdi_mobile + wdi_gerp + une_surlgpef + wdi_fertility, data=qog_std_ts_jan22, na.action=na.omit)

dem.rf
```

```
## 
## Call:
##  randomForest(formula = vdem_libdem ~ wdi_gendeqr + bci_bci +      wdi_gdpcapcon2010 + wdi_mobile + wdi_gerp + une_surlgpef +      wdi_fertility, data = qog_std_ts_jan22, na.action = na.omit) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 0.008876135
##                     % Var explained: 69.18
```

---

``` r
library(randomForestExplainer)
```

```
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
```

---

``` r
dem.mindepth <- min_depth_distribution(dem.rf)
plot_min_depth_distribution(dem.mindepth)
```

---

``` r
plot_min_depth_interactions(dem.rf)
```

```
## [1] "Warning: your forest does not contain information on local importance so 'mse_increase' measure cannot be extracted. To add it regrow the forest with the option localImp = TRUE and run this function again."
```

```
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## ℹ The deprecated feature was likely used in the randomForestExplainer package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
```

```
## Warning: There were 5600 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `bci_bci = min(bci_bci, na.rm = TRUE)`.
## ℹ In group 2: `tree = 1` and `split var = "une_surlgpef"`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 5599 remaining warnings.
```

---
### Causal Forests

Run a random forest predicting the treatment effect at each node, rather than the average value of the outcome.

Specifically, split groups up to maximize:

`$$n_{L} n_{R} (\hat{\tau} L − \hat{\tau} R)^2$$`

---
Robins, Rotnitzky & Zhao (1994) proved that the asymptotically optimal estimator for `$\tau$` is the Augmented Inverse Probability Weighted (AIPW) estimator.

---
<img src="images/AIPW.png" width="90%" style="display: block; margin: auto;" />

---

Using the AIPW to combine estimates that derive from causal forests gives us a *doubly robust* estimator.

---

```
## Warning: package 'grf' was built under R version 4.4.3
```

---

``` r
demdata <- qog_std_ts_jan22

demdata <- demdata %>% filter(!is.na(vdem_libdem) & !is.na(br_pres))

demdata.X <- with(demdata, rbind(wdi_gendeqr,bci_bci, wdi_gdpcapcon2010, wdi_mobile, wdi_gerp, une_surlgpef, wdi_fertility))

dem.cf <- causal_forest(Y=demdata$vdem_libdem, W=demdata$br_pres, X=t(demdata.X))

dem.cf
```

```
## GRF forest object of type causal_forest 
## Number of trees: 2000 
## Number of training samples: 9970 
## Variable importance: 
##     1     2     3     4     5     6     7 
## 0.009 0.148 0.330 0.279 0.017 0.020 0.197
```

---

``` r
demscores <- data.frame(score=get_scores(dem.cf))

p <- demscores %>%
  ggplot( aes(x=score)) +
  geom_histogram( binwidth=0.1, fill="#69b3a2", color="#e9ecef", alpha=0.9) +
  theme_ipsum() +
  theme(
    plot.title = element_text(size=15)
  )
```
---

``` r
p
```

```
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_bin()`).
```

```
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database
```

```
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
```

---
<img src="images/cfcausalscores.png" width="90%" style="display: block; margin: auto;" />

---

``` r
average_treatment_effect(dem.cf, target.sample = "overlap")
```

```
##     estimate      std.err 
## -0.046837070  0.005012132
```

``` r
average_treatment_effect(dem.cf, target.sample = "control")
```

```
##     estimate      std.err 
## -0.017043589  0.004399897
```

``` r
average_treatment_effect(dem.cf, target.sample = "treated")
```

```
## Warning in average_treatment_effect(dem.cf, target.sample = "treated"):
## Estimated treatment propensities go as high as 1 which means that treatment
## effects for some treated units may not be well identified. In this case, using
## `target.sample=control` may be helpful.
```

```
##     estimate      std.err 
## -0.041868636  0.004103802
```

---

``` r
temp.df <- data.frame(demscore = demscores$score, ht_colonial = as.factor(demdata$ht_colonial))

htcolonial_text <- c("Not colonized", "Dutch", "Spanish", "Italian", "US", "British", "French", "Portuguese", "Belgian", "Brit/French", "Austral.")

p2 <- temp.df %>% ggplot( aes(x=ht_colonial, y=demscore)) + 
    geom_violin() +
    xlab("Colonial History") +
    theme(legend.position="none") +
    xlab("") + scale_x_discrete(labels= htcolonial_text)
```

---

``` r
p2
```

```
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
```

---

``` r
temp.df <- data.frame(demscore = demscores$score, pop65 = demdata$wdi_pop65)

p3 <- temp.df %>% ggplot(aes(x=pop65, y=demscore)) +
  geom_point() +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
  theme_ipsum()
```

---

``` r
p3
```

```
## `geom_smooth()` using formula = 'y ~ x'
```

```
## Warning: Removed 914 rows containing non-finite outside the scale range
## (`stat_smooth()`).
```

```
## Warning: Removed 914 rows containing missing values or values outside the scale range
## (`geom_point()`).
```

```
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
```

---
<img src="images/cfpop65scatterplot.png" width="90%" style="display: block; margin: auto;" />

---
### Picking Control Variables

Throughout this course, we've often assumed we could correctly identify control variables. Can we, though?

---
### Picking Control Variables

Suppose we know enough to rule out post-treatment variables and colliders, but possibly not enough to distinguish a priori between confounders, irrelevant controls, and instruments. Suppose also that the total number of potential controls is very large. What to do?

---
### LASSO

Estimate a regression, as usual (minimizing the sum of squared errors), but subject to the added penalty term:

`$$\lambda \sum_{j=1}^{p}|\beta_{j}|$$`

---
### Picking Control Variables

If you are pretty sure you don't have any potential instruments in the data, and that there are a reasonable number of causes of the treatment that are not included in the dataset, then there is a double selection LASSO approach to picking variables.

---
<img src="images/Belloni1.png" width="90%" style="display: block; margin: auto;" />

---
### Double-selection

1. Run a LASSO regression predicting the outcome based on all the possible confounders.

2. Run a second LASSO regression predicting the treatment based on all the possible confounders.

3. Select all variables that have non-zero coefficients in either regression and use them as controls in your actual causal inference.

---

``` r
qogwdi <- qog_std_ts_jan22 %>% filter(year > 1979) %>% select(starts_with("wdi") | starts_with("year") | starts_with("cname"))
qogwdi <- qogwdi[,colSums(is.na(qogwdi))<3000]
qogwdi$vdem_libdem <- qog_std_ts_jan22$vdem_libdem[qog_std_ts_jan22$cname %in% qogwdi$cname & qog_std_ts_jan22$year >1979]
qogwdi.imp <- mice(qogwdi, method = "cart")
```

```
## 
##  iter imp variable
##   1   1  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   1   2  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   1   3  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   1   4  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   1   5  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   2   1  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   2   2  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   2   3  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   2   4  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   2   5  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   3   1  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   3   2  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   3   3  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   3   4  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   3   5  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   4   1  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   4   2  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   4   3  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   4   4  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   4   5  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   5   1  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   5   2  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   5   3  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   5   4  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
##   5   5  wdi_agedr  wdi_araland  wdi_area  wdi_birth  wdi_co2  wdi_death  wdi_enerenew  wdi_expmil  wdi_export  wdi_fdiin  wdi_fdiout  wdi_fertility  wdi_forest  wdi_gdpagr  wdi_gdpcapcon2010  wdi_gdpcapcur  wdi_gdpcapgr  wdi_gdpcappppcon2017  wdi_gdpcappppcur  wdi_gdpgr  wdi_gdpind  wdi_gdppppcon2017  wdi_gdppppcur  wdi_gerp  wdi_gerpf  wdi_gerpm  wdi_gnicapatlcur  wdi_gnicur  wdi_import  wdi_inflation  wdi_internet  wdi_lfpilo15  wdi_lifexp  wdi_lifexpf  wdi_lifexpm  wdi_mobile  wdi_mortf  wdi_mortinf  wdi_mortm  wdi_mortnn  wdi_mortu5  wdi_oilrent  wdi_pop  wdi_pop14  wdi_pop1564  wdi_pop65  wdi_popden  wdi_popf  wdi_popgr  wdi_poprul  wdi_poprulgr  wdi_popurbagr  wdi_refasy  wdi_refori  wdi_sva2010  wdi_svapg  wdi_svapgdp  wdi_tele  wdi_trade  wdi_tradeserv  wdi_wombuslawi  vdem_libdem
```

```
## Warning: Number of logged events: 1384
```
---

``` r
qogimputed1 <- mice::complete(qogwdi.imp)

summary(lm(wdi_gdpcapgr~vdem_libdem + I(log(wdi_gdpcapcon2010)) + wdi_gerp, data=qogimputed1))
```

```
## 
## Call:
## lm(formula = wdi_gdpcapgr ~ vdem_libdem + I(log(wdi_gdpcapcon2010)) + 
##     wdi_gerp, data = qogimputed1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.630  -2.211   0.260   2.530 138.901 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                0.747878   0.514770   1.453 0.146306    
## vdem_libdem                0.644974   0.302010   2.136 0.032740 *  
## I(log(wdi_gdpcapcon2010)) -0.079341   0.055517  -1.429 0.153002    
## wdi_gerp                   0.013138   0.003848   3.414 0.000644 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.279 on 8333 degrees of freedom
## Multiple R-squared:  0.002033,	Adjusted R-squared:  0.001674 
## F-statistic: 5.659 on 3 and 8333 DF,  p-value: 0.0007199
```

---

``` r
library(hdm)
```

---

``` r
demgrow.dualselection <- rlassoEffects(wdi_gdpcapgr ~ . - year - cname - cname_qog - cname_year - wdi_gdpgr - wdi_popgr, data=qogimputed1, I=~vdem_libdem)
demgrow.dualselection
```

```
## 
## Call:
## rlassoEffects.formula(formula = wdi_gdpcapgr ~ . - year - cname - 
##     cname_qog - cname_year - wdi_gdpgr - wdi_popgr, data = qogimputed1, 
##     I = ~vdem_libdem)
## 
## Coefficients:
## vdem_libdem  
##      0.1886
```

``` r
demgrow.dualselection$selection.matrix
```

```
##                      vdem_libdem
## wdi_agedr                  FALSE
## wdi_araland                FALSE
## wdi_area                   FALSE
## wdi_birth                   TRUE
## wdi_co2                    FALSE
## wdi_death                  FALSE
## wdi_enerenew               FALSE
## wdi_expmil                  TRUE
## wdi_export                  TRUE
## wdi_fdiin                  FALSE
## wdi_fdiout                 FALSE
## wdi_fertility               TRUE
## wdi_forest                  TRUE
## wdi_gdpagr                  TRUE
## wdi_gdpcapcon2010           TRUE
## wdi_gdpcapcur              FALSE
## wdi_gdpcappppcon2017       FALSE
## wdi_gdpcappppcur           FALSE
## wdi_gdpind                 FALSE
## wdi_gdppppcon2017          FALSE
## wdi_gdppppcur              FALSE
## wdi_gerp                   FALSE
## wdi_gerpf                  FALSE
## wdi_gerpm                  FALSE
## wdi_gniatlcur              FALSE
## wdi_gnicapatlcur            TRUE
## wdi_gnicur                 FALSE
## wdi_import                 FALSE
## wdi_inflation               TRUE
## wdi_internet                TRUE
## wdi_lfpilo15               FALSE
## wdi_lifexp                 FALSE
## wdi_lifexpf                FALSE
## wdi_lifexpm                FALSE
## wdi_mobile                 FALSE
## wdi_mortf                  FALSE
## wdi_mortinf                FALSE
## wdi_mortinff               FALSE
## wdi_mortinfm                TRUE
## wdi_mortm                  FALSE
## wdi_mortnn                 FALSE
## wdi_mortu5                 FALSE
## wdi_mortu5f                FALSE
## wdi_mortu5m                FALSE
## wdi_oilrent                 TRUE
## wdi_pop                    FALSE
## wdi_pop14                  FALSE
## wdi_pop1564                FALSE
## wdi_pop65                   TRUE
## wdi_popden                  TRUE
## wdi_popf                    TRUE
## wdi_poprul                  TRUE
## wdi_poprulgr               FALSE
## wdi_popurb                 FALSE
## wdi_popurbagr              FALSE
## wdi_refasy                 FALSE
## wdi_refori                 FALSE
## wdi_sva2010                FALSE
## wdi_svapg                   TRUE
## wdi_svapgdp                 TRUE
## wdi_tele                    TRUE
## wdi_trade                   TRUE
## wdi_tradeserv               TRUE
## wdi_wombuslawi              TRUE
## vdem_libdem                   NA
```

---
### Double-selection

One problem with double selection is that you will get instruments as control variables if they're in the data. Another problem is that it's inefficient because you aren't using information about the two relationships of interest jointly.

---
<img src="images/Shortreed1.png" width="90%" style="display: block; margin: auto;" />

---
### Adaptive LASSO

Adaptive LASSO allows the penalty term to change across coefficients:

`$$\lambda \sum_{j=1}^{p}|\hat{w}_{j}\beta_{j}|$$`

In a standard adaptive LASSO, `$\hat{w} = |\tilde{\beta}_{j}|^{-\gamma}$`, where `$\tilde{\beta}_{j}$` is a typical parametric estimate of the regression slopes and `$\gamma$` is a positive penalty parameter.

---
### Outcome-Adaptive LASSO

1. Run some kind of regression, predicting the outcome of interest with the treatment and all the potential control variables.

2. Run an adaptive lasso predicting the treatment based on the control variables. Set the penalty terms based on the coefficient estimates from step 1.

3. Use the control variables selected in step 2 to estimate a propensity score (or a regression or whatever).

---

``` r
library(glmnet)
```

```
## Loading required package: Matrix
```

```
## Loaded glmnet 4.1-8
```

---

``` r
qogimpcomplete <- qogimputed1[,colSums(is.na(qogimputed1))<1]

demgrow.adaptivelasso.part1 <- lm(wdi_gdpcapgr ~ . - year - cname - cname_qog - cname_year - wdi_gdpgr - wdi_popgr, data=qogimpcomplete)
tildebeta <- demgrow.adaptivelasso.part1$coefficients[2:59]
tildebeta[is.na(tildebeta)] <- .0000000000000000000000000000000000000000000000000000001
gammaval <- 0.5

xmat <- qogimpcomplete[,-c(17,20,49,62:66)]

demgrow.adaptivelasso.part2 <- glmnet(y=qogimputed1$wdi_gdpcapgr, x=xmat, penalty.factor=abs(tildebeta)^(-1*gammaval))
```

---

``` r
coef(demgrow.adaptivelasso.part2, 1)
```

```
## 59 x 1 sparse Matrix of class "dgCMatrix"
##                                 s1
## (Intercept)          -9.345746e+00
## wdi_agedr             5.534674e-02
## wdi_araland           1.373647e-03
## wdi_area              .           
## wdi_birth            -1.032587e-01
## wdi_co2              -3.521739e-02
## wdi_death             9.029891e-02
## wdi_enerenew          9.011291e-03
## wdi_expmil            1.410503e-02
## wdi_export            1.253821e-02
## wdi_fdiin             2.122848e-03
## wdi_fdiout            1.037552e-02
## wdi_fertility         2.376557e-01
## wdi_forest            4.922960e-03
## wdi_gdpagr           -3.871891e-02
## wdi_gdpcapcon2010     .           
## wdi_gdpcapcur        -4.365055e-06
## wdi_gdpcappppcon2017  .           
## wdi_gdpcappppcur     -8.084763e-06
## wdi_gdpind           -1.605287e-02
## wdi_gdppppcon2017     .           
## wdi_gdppppcur         .           
## wdi_gerp              5.182733e-02
## wdi_gerpf            -2.619747e-02
## wdi_gerpm            -2.373938e-02
## wdi_gnicapatlcur      .           
## wdi_gnicur            .           
## wdi_import            .           
## wdi_inflation        -9.724584e-04
## wdi_internet         -8.993308e-03
## wdi_lfpilo15         -6.833171e-03
## wdi_lifexp           -8.590309e-02
## wdi_lifexpf           .           
## wdi_lifexpm           8.262990e-02
## wdi_mobile           -1.648997e-03
## wdi_mortf            -2.125031e-03
## wdi_mortinf          -3.006591e-02
## wdi_mortm            -1.016983e-03
## wdi_mortnn            3.832627e-02
## wdi_mortu5            5.844230e-07
## wdi_oilrent           3.047363e-02
## wdi_pop               .           
## wdi_pop14            -7.895719e-02
## wdi_pop1564           9.824204e-02
## wdi_pop65            -1.526973e-01
## wdi_popden            1.517555e-04
## wdi_popf              1.709450e-01
## wdi_poprul            1.470778e-02
## wdi_poprulgr         -5.620704e-03
## wdi_popurbagr        -1.264933e-01
## wdi_refasy            .           
## wdi_refori            .           
## wdi_sva2010           .           
## wdi_svapg             3.266980e-01
## wdi_svapgdp          -4.406059e-02
## wdi_tele              .           
## wdi_trade            -7.087592e-03
## wdi_tradeserv         9.052865e-03
## wdi_wombuslawi        5.288284e-03
```
---

``` r
coef(demgrow.adaptivelasso.part2, 2)
```

``` r
coef(demgrow.adaptivelasso.part2, 5)
```

``` r
coef(demgrow.adaptivelasso.part2, 20)
```

``` r
coef(demgrow.adaptivelasso.part2, 100)
```

```
## 59 x 1 sparse Matrix of class "dgCMatrix"
##                                 s1
## (Intercept)          -7.8852698479
## wdi_agedr             0.0380252922
## wdi_araland           0.0010582689
## wdi_area              .           
## wdi_birth            -0.1008766799
## wdi_co2              -0.0471022393
## wdi_death             0.0939185445
## wdi_enerenew          0.0079707443
## wdi_expmil            0.0126511418
## wdi_export            0.0074040488
## wdi_fdiin             0.0018750691
## wdi_fdiout            0.0090568283
## wdi_fertility         0.2303643697
## wdi_forest            0.0045958984
## wdi_gdpagr           -0.0356760398
## wdi_gdpcapcon2010     .           
## wdi_gdpcapcur         .           
## wdi_gdpcappppcon2017  .           
## wdi_gdpcappppcur      .           
## wdi_gdpind           -0.0151387720
## wdi_gdppppcon2017     .           
## wdi_gdppppcur         .           
## wdi_gerp              0.0211437388
## wdi_gerpf            -0.0112513401
## wdi_gerpm            -0.0080292673
## wdi_gnicapatlcur      .           
## wdi_gnicur            .           
## wdi_import           -0.0021360434
## wdi_inflation        -0.0009520943
## wdi_internet         -0.0118452781
## wdi_lfpilo15         -0.0060847740
## wdi_lifexp           -0.0757632329
## wdi_lifexpf           .           
## wdi_lifexpm           0.0786970197
## wdi_mobile           -0.0008612862
## wdi_mortf            -0.0020418501
## wdi_mortinf          -0.0300044604
## wdi_mortm            -0.0005842015
## wdi_mortnn            0.0380866007
## wdi_mortu5            .           
## wdi_oilrent           0.0279754542
## wdi_pop               .           
## wdi_pop14            -0.0616318556
## wdi_pop1564           0.0686836984
## wdi_pop65            -0.1423627426
## wdi_popden            0.0001091097
## wdi_popf              0.1753770788
## wdi_poprul            0.0147829331
## wdi_poprulgr         -0.0056146201
## wdi_popurbagr        -0.1294318655
## wdi_refasy            .           
## wdi_refori            .           
## wdi_sva2010           .           
## wdi_svapg             0.3265806570
## wdi_svapgdp          -0.0435683960
## wdi_tele              .           
## wdi_trade            -0.0032744835
## wdi_tradeserv         0.0080185486
## wdi_wombuslawi        0.0036530968
```

---

``` r
qogimpcomplete <- qogimputed1[,colSums(is.na(qogimputed1))<1]

xmat <- qogimpcomplete[,-c(17,20,49,62:66)]

demgrow.adaptivelasso.part2 <- glmnet(y=qogimputed1$wdi_gdpcapgr, x=xmat, penalty.factor=abs(tildebeta)^(-1*gammaval))
```

---

``` r
coef(demgrow.adaptivelasso.part2, 1)
```

```
## 59 x 1 sparse Matrix of class "dgCMatrix"
##                                s1
## (Intercept)          -8.268600355
## wdi_agedr             .          
## wdi_araland           .          
## wdi_area              .          
## wdi_birth            -0.100726910
## wdi_co2              -0.030405643
## wdi_death             0.054934931
## wdi_enerenew          .          
## wdi_expmil            .          
## wdi_export            .          
## wdi_fdiin             .          
## wdi_fdiout            .          
## wdi_fertility         0.213900619
## wdi_forest            .          
## wdi_gdpagr           -0.002673310
## wdi_gdpcapcon2010     .          
## wdi_gdpcapcur         .          
## wdi_gdpcappppcon2017  .          
## wdi_gdpcappppcur      .          
## wdi_gdpind            .          
## wdi_gdppppcon2017     .          
## wdi_gdppppcur         .          
## wdi_gerp              0.002086477
## wdi_gerpf             .          
## wdi_gerpm             .          
## wdi_gnicapatlcur      .          
## wdi_gnicur            .          
## wdi_import            .          
## wdi_inflation         .          
## wdi_internet          .          
## wdi_lfpilo15          .          
## wdi_lifexp           -0.029147464
## wdi_lifexpf           .          
## wdi_lifexpm           0.061713464
## wdi_mobile            .          
## wdi_mortf             .          
## wdi_mortinf          -0.001893674
## wdi_mortm             .          
## wdi_mortnn            .          
## wdi_mortu5            .          
## wdi_oilrent           0.008902340
## wdi_pop               .          
## wdi_pop14            -0.001117606
## wdi_pop1564           0.005267741
## wdi_pop65            -0.101444663
## wdi_popden            .          
## wdi_popf              0.195803214
## wdi_poprul            .          
## wdi_poprulgr          .          
## wdi_popurbagr        -0.132569663
## wdi_refasy            .          
## wdi_refori            .          
## wdi_sva2010           .          
## wdi_svapg             0.332691655
## wdi_svapgdp          -0.022511383
## wdi_tele              .          
## wdi_trade             .          
## wdi_tradeserv         .          
## wdi_wombuslawi        .
```
---

``` r
coef(demgrow.adaptivelasso.part2, 2)
```

``` r
coef(demgrow.adaptivelasso.part2, 5)
```

``` r
coef(demgrow.adaptivelasso.part2, 20)
```

``` r
coef(demgrow.adaptivelasso.part2, 100)
```

---
### Textual Data and Causal Inference