R for Data Analysis: A Short Tutorial

class: top, left, inverse, title-slide

.title[
# R for Data Analysis: A Short Tutorial
]
.subtitle[
## Session 3: Modelling data
]
.author[
### Dimiter Toshkov
]
.institute[
### Institute of Public Administration, Leiden University
]
.date[
### last updated: 2025-03-31
]

---

# Last session...

.title-slide {
  background-image: url(https://cran.r-project.org/Rlogo.svg);
  background-position: 50% 0%; ## just start changing this
  background-size: 150px;
  background-color: #fff;
  padding-left: 100px;  /* delete this for 4:3 aspect ratio */
}

.remark-slide-content {
    font-size: 28px;
    padding: 1em 1em 1em 1em;
}

.remark-slide-content > h1 {
  font-size: 32px;
  margin-top: -85px;
}

</style>

---
# Ready for some modelling?

.center[![trachoma](data:image/png;base64,#figs/zoolander2.jpg)]

---
# Let's load some data!

``` r
library(tidyverse)
library(haven)
dt <- read_sav("./data/ESS11.sav")
head(dt)
## # A tibble: 6 × 640
##   name     essround edition proddate    idno cntry       dweight pspwght pweight
##   <chr>       <dbl> <chr>   <chr>      <dbl> <chr+lbl>     <dbl>   <dbl>   <dbl>
## 1 ESS11e02       11 2.0     20.11.2024 50014 AT [Austri…   1.19    0.393   0.331
## 2 ESS11e02       11 2.0     20.11.2024 50030 AT [Austri…   0.610   0.325   0.331
## 3 ESS11e02       11 2.0     20.11.2024 50057 AT [Austri…   1.39    4.00    0.331
## 4 ESS11e02       11 2.0     20.11.2024 50106 AT [Austri…   0.556   0.176   0.331
## 5 ESS11e02       11 2.0     20.11.2024 50145 AT [Austri…   0.723   1.06    0.331
## 6 ESS11e02       11 2.0     20.11.2024 50158 AT [Austri…   0.993   1.39    0.331
## # ℹ 631 more variables: anweight <dbl>, nwspol <dbl+lbl>, netusoft <dbl+lbl>,
## #   netustm <dbl+lbl>, ppltrst <dbl+lbl>, pplfair <dbl+lbl>, pplhlp <dbl+lbl>,
## #   polintr <dbl+lbl>, psppsgva <dbl+lbl>, actrolga <dbl+lbl>,
## #   psppipla <dbl+lbl>, cptppola <dbl+lbl>, trstprl <dbl+lbl>,
## #   trstlgl <dbl+lbl>, trstplc <dbl+lbl>, trstplt <dbl+lbl>, trstprt <dbl+lbl>,
## #   trstep <dbl+lbl>, trstun <dbl+lbl>, vote <dbl+lbl>, prtvtdat <dbl+lbl>,
## #   prtvtebe <dbl+lbl>, prtvtchr <dbl+lbl>, prtvtccy <dbl+lbl>, …
```

---
# Plotting distributions (1)
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M396.8 352h22.4c6.4 0 12.8-6.4 12.8-12.8V108.8c0-6.4-6.4-12.8-12.8-12.8h-22.4c-6.4 0-12.8 6.4-12.8 12.8v230.4c0 6.4 6.4 12.8 12.8 12.8zm-192 0h22.4c6.4 0 12.8-6.4 12.8-12.8V140.8c0-6.4-6.4-12.8-12.8-12.8h-22.4c-6.4 0-12.8 6.4-12.8 12.8v198.4c0 6.4 6.4 12.8 12.8 12.8zm96 0h22.4c6.4 0 12.8-6.4 12.8-12.8V204.8c0-6.4-6.4-12.8-12.8-12.8h-22.4c-6.4 0-12.8 6.4-12.8 12.8v134.4c0 6.4 6.4 12.8 12.8 12.8zM496 400H48V80c0-8.84-7.16-16-16-16H16C7.16 64 0 71.16 0 80v336c0 17.67 14.33 32 32 32h464c8.84 0 16-7.16 16-16v-16c0-8.84-7.16-16-16-16zm-387.2-48h22.4c6.4 0 12.8-6.4 12.8-12.8v-70.4c0-6.4-6.4-12.8-12.8-12.8h-22.4c-6.4 0-12.8 6.4-12.8 12.8v70.4c0 6.4 6.4 12.8 12.8 12.8z"></path></svg> Before we decide on a statistical model, we want to inspect the distribution of the variables, for example by looking at histograms:

.pull-left[

``` r
#plot first variable
hist(dt$trstep, 
     col=rgb(0,0,1,0.25), 
     breaks=10,
     xlab='Trust in the European Parliament (blue)\n
     and in national parliaments (red)',
     main=''
) 
#add second variable
hist(dt$trstprl, 
     col=rgb(1,0,0,0.25), 
     breaks=10,
     add=T) 
```
]

.pull-right[
<img src="data:image/png;base64,#R-Tutorial-Session-3_files/figure-html/unnamed-chunk-3-1.png" width="504" />
]

---
# Plotting distributions (2) 
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M500 384c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12H12c-6.6 0-12-5.4-12-12V76c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v308h436zM372.7 159.5L288 216l-85.3-113.7c-5.1-6.8-15.5-6.3-19.9 1L96 248v104h384l-89.9-187.8c-3.2-6.5-11.4-8.7-17.4-4.7z"></path></svg> Histograms show frequencies. If we want to see relative proportions, we switch to densities 
.pull-left[

``` r
# plot first variable
plot(density(
  dt$agea[dt$gndr == 1], na.rm=TRUE), 
  col=rgb(0, 0, 1, 0.75), # semi-transparent colors
  xlab = 'Age (in years)',
  main = "") 
# add second
lines(density(
  dt$agea[dt$gndr == 2], na.rm=TRUE), 
   col = rgb(1, 0, 0, 0.75)
  ) 
```
]
.pull-right[
<img src="data:image/png;base64,#R-Tutorial-Session-3_files/figure-html/unnamed-chunk-5-1.png" width="504" />
]

---
# Cross-tabulations
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M464 32H48C21.49 32 0 53.49 0 80v352c0 26.51 21.49 48 48 48h416c26.51 0 48-21.49 48-48V80c0-26.51-21.49-48-48-48zM224 416H64v-96h160v96zm0-160H64v-96h160v96zm224 160H288v-96h160v96zm0-160H288v-96h160v96z"></path></svg> Cross-tabs are a common first step to examine associations between (categorical) variables.

``` r
dt$voted.dummy = ifelse(dt$vote == 1, 1, 0)
dt$sex = ifelse (dt$gndr == 1, 'male', 'female')
table (dt$voted.dummy, dt$sex) # frequencies
##    
##     female  male
##   0   5581  4664
##   1  15607 13944
table (dt$voted.dummy, dt$sex) %>% # percentages
  prop.table(2) %>% round(2)
##    
##     female male
##   0   0.26 0.25
##   1   0.74 0.75
```

---
# Tables of descriptive stats
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M464 32H48C21.49 32 0 53.49 0 80v352c0 26.51 21.49 48 48 48h416c26.51 0 48-21.49 48-48V80c0-26.51-21.49-48-48-48zM224 416H64v-96h160v96zm0-160H64v-96h160v96zm224 160H288v-96h160v96zm0-160H288v-96h160v96z"></path></svg> Often, before we report any results from statistical models, we want to present a descriptive summary of the variables used. We already know the `summary()` function, which we can use to do exactly that, but its output is not very easy to insert into a text file. 
The `modelsummary` package has some great functions that produce tables we can easily insert in `MSWord` and other types of files. For details, see [here](https://vincentarelbundock.github.io/modelsummary/articles/datasummary.html). Another great package for tables of descriptives is `gtsummary`.

---
# Table of descriptive stats: tbl_summary()

``` r
library(gtsummary)
dt %>% dplyr::select(stfdem, lrscale, netustm, sex) %>% tbl_summary(by='sex')
```

<div id="ihkqtrcfsr" style="padding-left:0px;padding-right:0px;padding-top:10px;padding-bottom:10px;overflow-x:auto;overflow-y:auto;width:auto;height:auto;">
<style>#ihkqtrcfsr table {
  font-family: system-ui, 'Segoe UI', Roboto, Helvetica, Arial, sans-serif, 'Apple Color Emoji', 'Segoe UI Emoji', 'Segoe UI Symbol', 'Noto Color Emoji';
  -webkit-font-smoothing: antialiased;
  -moz-osx-font-smoothing: grayscale;
}

#ihkqtrcfsr thead, #ihkqtrcfsr tbody, #ihkqtrcfsr tfoot, #ihkqtrcfsr tr, #ihkqtrcfsr td, #ihkqtrcfsr th {
  border-style: none;
}

#ihkqtrcfsr p {
  margin: 0;
  padding: 0;
}

#ihkqtrcfsr .gt_table {
  display: table;
  border-collapse: collapse;
  line-height: normal;
  margin-left: auto;
  margin-right: auto;
  color: #333333;
  font-size: 16px;
  font-weight: normal;
  font-style: normal;
  background-color: #FFFFFF;
  width: auto;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #A8A8A8;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #A8A8A8;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
}

#ihkqtrcfsr .gt_caption {
  padding-top: 4px;
  padding-bottom: 4px;
}

#ihkqtrcfsr .gt_title {
  color: #333333;
  font-size: 125%;
  font-weight: initial;
  padding-top: 4px;
  padding-bottom: 4px;
  padding-left: 5px;
  padding-right: 5px;
  border-bottom-color: #FFFFFF;
  border-bottom-width: 0;
}

#ihkqtrcfsr .gt_subtitle {
  color: #333333;
  font-size: 85%;
  font-weight: initial;
  padding-top: 3px;
  padding-bottom: 5px;
  padding-left: 5px;
  padding-right: 5px;
  border-top-color: #FFFFFF;
  border-top-width: 0;
}

#ihkqtrcfsr .gt_heading {
  background-color: #FFFFFF;
  text-align: center;
  border-bottom-color: #FFFFFF;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
}

#ihkqtrcfsr .gt_bottom_border {
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#ihkqtrcfsr .gt_col_headings {
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
}

#ihkqtrcfsr .gt_col_heading {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: normal;
  text-transform: inherit;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: bottom;
  padding-top: 5px;
  padding-bottom: 6px;
  padding-left: 5px;
  padding-right: 5px;
  overflow-x: hidden;
}

#ihkqtrcfsr .gt_column_spanner_outer {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: normal;
  text-transform: inherit;
  padding-top: 0;
  padding-bottom: 0;
  padding-left: 4px;
  padding-right: 4px;
}

#ihkqtrcfsr .gt_column_spanner_outer:first-child {
  padding-left: 0;
}

#ihkqtrcfsr .gt_column_spanner_outer:last-child {
  padding-right: 0;
}

#ihkqtrcfsr .gt_column_spanner {
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  vertical-align: bottom;
  padding-top: 5px;
  padding-bottom: 5px;
  overflow-x: hidden;
  display: inline-block;
  width: 100%;
}

#ihkqtrcfsr .gt_spanner_row {
  border-bottom-style: hidden;
}

#ihkqtrcfsr .gt_group_heading {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: middle;
  text-align: left;
}

#ihkqtrcfsr .gt_empty_group_heading {
  padding: 0.5px;
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  vertical-align: middle;
}

#ihkqtrcfsr .gt_from_md > :first-child {
  margin-top: 0;
}

#ihkqtrcfsr .gt_from_md > :last-child {
  margin-bottom: 0;
}

#ihkqtrcfsr .gt_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  margin: 10px;
  border-top-style: solid;
  border-top-width: 1px;
  border-top-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: middle;
  overflow-x: hidden;
}

#ihkqtrcfsr .gt_stub {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-right-style: solid;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  padding-left: 5px;
  padding-right: 5px;
}

#ihkqtrcfsr .gt_stub_row_group {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-right-style: solid;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  padding-left: 5px;
  padding-right: 5px;
  vertical-align: top;
}

#ihkqtrcfsr .gt_row_group_first td {
  border-top-width: 2px;
}

#ihkqtrcfsr .gt_row_group_first th {
  border-top-width: 2px;
}

#ihkqtrcfsr .gt_summary_row {
  color: #333333;
  background-color: #FFFFFF;
  text-transform: inherit;
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
}

#ihkqtrcfsr .gt_first_summary_row {
  border-top-style: solid;
  border-top-color: #D3D3D3;
}

#ihkqtrcfsr .gt_first_summary_row.thick {
  border-top-width: 2px;
}

#ihkqtrcfsr .gt_last_summary_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#ihkqtrcfsr .gt_grand_summary_row {
  color: #333333;
  background-color: #FFFFFF;
  text-transform: inherit;
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
}

#ihkqtrcfsr .gt_first_grand_summary_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-top-style: double;
  border-top-width: 6px;
  border-top-color: #D3D3D3;
}

#ihkqtrcfsr .gt_last_grand_summary_row_top {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-bottom-style: double;
  border-bottom-width: 6px;
  border-bottom-color: #D3D3D3;
}

#ihkqtrcfsr .gt_striped {
  background-color: rgba(128, 128, 128, 0.05);
}

#ihkqtrcfsr .gt_table_body {
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#ihkqtrcfsr .gt_footnotes {
  color: #333333;
  background-color: #FFFFFF;
  border-bottom-style: none;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
}

#ihkqtrcfsr .gt_footnote {
  margin: 0px;
  font-size: 90%;
  padding-top: 4px;
  padding-bottom: 4px;
  padding-left: 5px;
  padding-right: 5px;
}

#ihkqtrcfsr .gt_sourcenotes {
  color: #333333;
  background-color: #FFFFFF;
  border-bottom-style: none;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
}

#ihkqtrcfsr .gt_sourcenote {
  font-size: 90%;
  padding-top: 4px;
  padding-bottom: 4px;
  padding-left: 5px;
  padding-right: 5px;
}

#ihkqtrcfsr .gt_left {
  text-align: left;
}

#ihkqtrcfsr .gt_center {
  text-align: center;
}

#ihkqtrcfsr .gt_right {
  text-align: right;
  font-variant-numeric: tabular-nums;
}

#ihkqtrcfsr .gt_font_normal {
  font-weight: normal;
}

#ihkqtrcfsr .gt_font_bold {
  font-weight: bold;
}

#ihkqtrcfsr .gt_font_italic {
  font-style: italic;
}

#ihkqtrcfsr .gt_super {
  font-size: 65%;
}

#ihkqtrcfsr .gt_footnote_marks {
  font-size: 75%;
  vertical-align: 0.4em;
  position: initial;
}

#ihkqtrcfsr .gt_asterisk {
  font-size: 100%;
  vertical-align: 0;
}

#ihkqtrcfsr .gt_indent_1 {
  text-indent: 5px;
}

#ihkqtrcfsr .gt_indent_2 {
  text-indent: 10px;
}

#ihkqtrcfsr .gt_indent_3 {
  text-indent: 15px;
}

#ihkqtrcfsr .gt_indent_4 {
  text-indent: 20px;
}

#ihkqtrcfsr .gt_indent_5 {
  text-indent: 25px;
}

#ihkqtrcfsr .katex-display {
  display: inline-flex !important;
  margin-bottom: 0.75em !important;
}

#ihkqtrcfsr div.Reactable > div.rt-table > div.rt-thead > div.rt-tr.rt-tr-group-header > div.rt-th-group:after {
  height: 0px !important;
}
</style>
<table class="gt_table" data-quarto-disable-processing="false" data-quarto-bootstrap="false">
  <thead>
    <tr class="gt_col_headings">
      <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1" scope="col" id="label"><span class='gt_from_md'><strong>Characteristic</strong></span></th>
      <th class="gt_col_heading gt_columns_bottom_border gt_center" rowspan="1" colspan="1" scope="col" id="stat_1"><span class='gt_from_md'><strong>female</strong><br />
N = 21,396</span><span class="gt_footnote_marks" style="white-space:nowrap;font-style:italic;font-weight:normal;line-height:0;"><sup>1</sup></span></th>
      <th class="gt_col_heading gt_columns_bottom_border gt_center" rowspan="1" colspan="1" scope="col" id="stat_2"><span class='gt_from_md'><strong>male</strong><br />
N = 18,760</span><span class="gt_footnote_marks" style="white-space:nowrap;font-style:italic;font-weight:normal;line-height:0;"><sup>1</sup></span></th>
    </tr>
  </thead>
  <tbody class="gt_table_body">
    <tr><td headers="label" class="gt_row gt_left">How satisfied with the way democracy works in country</td>
<td headers="stat_1" class="gt_row gt_center">5.00 (3.00, 7.00)</td>
<td headers="stat_2" class="gt_row gt_center">5.00 (3.00, 7.00)</td></tr>
    <tr><td headers="label" class="gt_row gt_left">    Unknown</td>
<td headers="stat_1" class="gt_row gt_center">846</td>
<td headers="stat_2" class="gt_row gt_center">398</td></tr>
    <tr><td headers="label" class="gt_row gt_left">Placement on left right scale</td>
<td headers="stat_1" class="gt_row gt_center">5.00 (4.00, 6.00)</td>
<td headers="stat_2" class="gt_row gt_center">5.00 (4.00, 7.00)</td></tr>
    <tr><td headers="label" class="gt_row gt_left">    Unknown</td>
<td headers="stat_1" class="gt_row gt_center">3,017</td>
<td headers="stat_2" class="gt_row gt_center">1,909</td></tr>
    <tr><td headers="label" class="gt_row gt_left">Internet use, how much time on typical day, in minutes</td>
<td headers="stat_1" class="gt_row gt_center">180 (90, 300)</td>
<td headers="stat_2" class="gt_row gt_center">180 (120, 330)</td></tr>
    <tr><td headers="label" class="gt_row gt_left">    Unknown</td>
<td headers="stat_1" class="gt_row gt_center">4,940</td>
<td headers="stat_2" class="gt_row gt_center">3,818</td></tr>
  </tbody>
  
  <tfoot class="gt_footnotes">
    <tr>
      <td class="gt_footnote" colspan="3"><span class="gt_footnote_marks" style="white-space:nowrap;font-style:italic;font-weight:normal;line-height:0;"><sup>1</sup></span> <span class='gt_from_md'>Median (Q1, Q3)</span></td>
    </tr>
  </tfoot>
</table>
</div>

---
# Table of descriptive stats: datasummary_skim()

``` r
library(modelsummary)
dt.s<-data.frame(select(dt, lrscale, trstun, trstep, trstprl))
datasummary_skim(dt.s, output = 'markdown')
# change the type of output to suit your desired destination
```

---
# Summarize per group
But cross-tabs aren't so great when we want to summarize a continuous variable per group. Then, we can then use `summarize()`. There are [alternatives](https://stackoverflow.com/questions/9847054/how-to-get-summary-statistics-by-group) as well.

``` r
dt %>% 
  group_by(sex) %>% 
  summarize(mean = mean(trstun, na.rm=TRUE),
            sd = sd(trstun, na.rm=TRUE),
            min = min(trstun, na.rm=TRUE),
            max = max(trstun, na.rm=TRUE))
## # A tibble: 2 × 5
##   sex     mean    sd min                 max                
##   <chr>  <dbl> <dbl> <dbl+lbl>           <dbl+lbl>          
## 1 female  5.18  2.57 0 [No trust at all] 10 [Complete trust]
## 2 male    5.06  2.62 0 [No trust at all] 10 [Complete trust]
```

---
# OK, OK, let's do a **t-test**

Two-sample t-test of *Trust in UN* per *Gender*:

``` r
t.test(
  dt$trstun[dt$gndr == 1], 
  dt$trstun[dt$gndr == 2])
## 
## 	Welch Two Sample t-test
## 
## data:  dt$trstun[dt$gndr == 1] and dt$trstun[dt$gndr == 2]
## t = -4.2864, df = 36983, p-value = 1.821e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.16781758 -0.06250018
## sample estimates:
## mean of x mean of y 
##  5.061234  5.176393
```

---
# Yes, we can also do chi-squared tests
Chi-squared (χ2) tests are appropriate when we want to examine associations between categorical (factor) variables and for contingency tables.

Voting per gender:

``` r
chisq.test(factor(dt$gndr), factor(dt$vote))
## 
## 	Pearson's Chi-squared test
## 
## data:  factor(dt$gndr) and factor(dt$vote)
## X-squared = 38.022, df = 2, p-value = 5.543e-09
```

---
# Correlations, correlations
We compute simple bivariate correlations with `cor(x,y, use='complete')`. But the `corrplot` package provides for nice correlation plots: [details](https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html).

.panelset[
.panel[.panel-name[Corrplot code]

``` r
# first, subset to drop missing data
dt.com <- dt %>% 
  select(lrscale, trstep, trstun, trstprl, stfdem,  stflife) %>% 
  drop_na()

# then, run corrplot
library(corrplot)
corrplot.mixed(cor(dt.com), lower.col = "black", number.cex = .7)
```
]
.panel[.panel-name[Resulting plot]
<img src="data:image/png;base64,#R-Tutorial-Session-3_files/figure-html/unnamed-chunk-13-1.png" width="504" />

]
]

---
# Good old OLS: Linear regression
Naturally, our first stat model should be the workhorse of data analysis, linear regression. Note that we assign the output to a named object and then use `summary()` on the object to get an overview of the results (scroll down the slide to see all results)

``` r
model.lm.1 <- lm(stfdem ~ lrscale + trstprl + factor(sex) + netustm + etfruit + cntry, 
               weights=pspwght,
               data=dt)
summary(model.lm.1)
```

```
## 
## Call:
## lm(formula = stfdem ~ lrscale + trstprl + factor(sex) + netustm + 
##     etfruit + cntry, data = dt, weights = pspwght)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5250  -1.1342   0.0926   1.1351  14.4580 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.603e+00  6.639e-02  39.210  < 2e-16 ***
## lrscale          8.058e-02  5.305e-03  15.190  < 2e-16 ***
## trstprl          5.547e-01  4.851e-03 114.335  < 2e-16 ***
## factor(sex)male  6.490e-02  2.338e-02   2.775 0.005516 ** 
## netustm         -7.609e-05  6.221e-05  -1.223 0.221301    
## etfruit         -3.300e-02  9.141e-03  -3.610 0.000307 ***
## cntryBE         -2.785e-01  7.128e-02  -3.907 9.35e-05 ***
## cntryCH          1.062e+00  7.449e-02  14.261  < 2e-16 ***
## cntryCY         -6.786e-01  1.035e-01  -6.554 5.70e-11 ***
## cntryDE          4.393e-02  6.389e-02   0.688 0.491718    
## cntryES         -5.586e-01  7.028e-02  -7.948 1.97e-15 ***
## cntryFI          2.741e-01  7.074e-02   3.875 0.000107 ***
## cntryFR         -9.848e-01  7.185e-02 -13.707  < 2e-16 ***
## cntryGB         -4.608e-01  7.275e-02  -6.335 2.41e-10 ***
## cntryGR         -4.858e-01  6.578e-02  -7.385 1.57e-13 ***
## cntryHR         -1.457e-01  8.002e-02  -1.821 0.068583 .  
## cntryHU         -1.148e+00  7.136e-02 -16.081  < 2e-16 ***
## cntryIE          2.654e-01  6.984e-02   3.800 0.000145 ***
## cntryIS         -1.547e-02  8.450e-02  -0.183 0.854780    
## cntryIT         -3.260e-01  6.851e-02  -4.758 1.97e-06 ***
## cntryLT         -2.802e-01  8.321e-02  -3.368 0.000759 ***
## cntryNL          1.426e-01  6.850e-02   2.081 0.037442 *  
## cntryNO          6.581e-01  7.283e-02   9.037  < 2e-16 ***
## cntryPL          1.455e-02  8.010e-02   0.182 0.855899    
## cntryPT         -6.793e-01  8.231e-02  -8.253  < 2e-16 ***
## cntryRS         -1.506e+00  8.507e-02 -17.708  < 2e-16 ***
## cntrySE          1.792e-01  7.480e-02   2.396 0.016583 *  
## cntrySI         -2.609e-01  8.322e-02  -3.135 0.001721 ** 
## cntrySK         -1.808e-01  8.054e-02  -2.245 0.024779 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.925 on 27501 degrees of freedom
##   (12626 observations deleted due to missingness)
## Multiple R-squared:  0.4324,	Adjusted R-squared:  0.4318 
## F-statistic: 748.1 on 28 and 27501 DF,  p-value: < 2.2e-16
```

---
# Running statistical models in R (1)
There is a common template behind the different stat models. We put the outcome variable on the left side, then we put tilde (`~`), then we put the predictors separated by `+`. After we are done, we put a comma (` ,`) and specify the `data` argument, and other optional arguments such as `weights` or distribution family. An intercept is usually included by default. We can include interactions with `*`.

We can store the result of the estimation in a (named) object, which would be of type `list`. Then we can call functions, such as `print()`, `summary()` and `plot()` on this object. 
---
# Running statistical models in R (2)
We can also access particular elements of this list. For example, for `lm` models, `model.lm.1[[1]]` will give the estimated coefficients, as well `model.lm.1$coef`.

We can access the residuals with `model.lm.1$residuals` and predicted values with `predict(model.lm.1)` (we can also supply new values for the predictors for which to compute predictions).

---
# Diagnostics
.pull-left[
If you are into these things, you can get some pre-fabricated diagnostics for your linear models with `plot(model.lm.1)` or with the package `gglm` and `gglm()`.]

.pull-right[
<img src="data:image/png;base64,#R-Tutorial-Session-3_files/figure-html/unnamed-chunk-15-1.png" width="504" />
]

---
# What if you care about being robust?

You can make your linear regression provide robust or clustered standard errors. This used to be a hassle, and a major reason for people to stick to other proprietary data-analytic software. But not anymore. Use the `fixest` [package](https://cran.r-project.org/web/packages/fixest/vignettes/fixest_walkthrough.html) . The syntax is very similar to `lm`. But we specify fixed effects slightly differently by using the `|` symbol.

Alternatively, we can use the `lm_robust()` function from the `estimatr` [package](https://declaredesign.org/r/estimatr/reference/lm_robust.html).

Andrew Heiss has [a great overview](https://evalf21.classes.andrewheiss.com/example/standard-errors/) of the options.

---
# Robust errors in action
An example of a linear regression with fixed effects at the country level and heteroskedastic robust standard errors with the `feols()` function from `fixest`. Using this model we make the amazing discovery that eating more fruit (`etfruit`) decreases *significantly* satisfaction with democracy! **(Scroll down the slide for full results.)**

``` r
library(fixest)
model.lm.1.robust<-feols(stfdem ~ lrscale + trstprl + factor(sex) + 
                           netustm + etfruit | cntry, # fixed effects per country
                         vcov = "hetero", # heteroskedastic robust SEs
                         data = dt)
summary(model.lm.1.robust)
## OLS estimation, Dep. Var.: stfdem
## Observations: 27,530
## Fixed-effects: cntry: 24
## Standard-errors: Heteroskedasticity-robust 
##                  Estimate Std. Error   t value   Pr(>|t|)    
## lrscale          0.080249   0.005882  13.64397  < 2.2e-16 ***
## trstprl          0.550878   0.005501 100.14107  < 2.2e-16 ***
## factor(sex)male  0.056556   0.023450   2.41175 1.5883e-02 *  
## netustm         -0.000067   0.000063  -1.06095 2.8872e-01    
## etfruit         -0.040929   0.009549  -4.28628 1.8231e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.91071     Adj. R2: 0.427321
##                 Within R2: 0.326625
```

---
# Generalized linear models
The `glm()` function has you covered for various types of models: 
- logistic (`glm ( y ~ x, family = 'binomial')`) for binary data, 
- Poisson (`glm ( y ~ x, family = 'poisson')`) for count data, 
- negative binomial for overdispersed count data (`glm( y ~ x, family = 'quasipoisson')`), and more.

---
# Logistic regression results
Let's see the `glm()` function in action performing a simple logistic regression:

``` r
model.glm.1<-glm(voted.dummy ~ stfdem + sex + lrscale + trstprl, 
                 data = dt, family = 'binomial')
exp(model.glm.1$coefficients) # to exponentiate the coefficients
## (Intercept)      stfdem     sexmale     lrscale     trstprl 
##   2.5651376   0.9844047   1.0179269   1.0196026   1.0693732
summary(model.glm.1)
## 
## Call:
## glm(formula = voted.dummy ~ stfdem + sex + lrscale + trstprl, 
##     family = "binomial", data = dt)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9633   0.6193   0.6762   0.7215   0.8613  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.942012   0.041820  22.525  < 2e-16 ***
## stfdem      -0.015718   0.006568  -2.393 0.016702 *  
## sexmale      0.017768   0.026318   0.675 0.499597    
## lrscale      0.019413   0.005895   3.293 0.000991 ***
## trstprl      0.067073   0.006343  10.574  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 35961  on 34270  degrees of freedom
## Residual deviance: 35808  on 34266  degrees of freedom
##   (5885 observations deleted due to missingness)
## AIC: 35818
## 
## Number of Fisher Scoring iterations: 4
```

---
# Models with interactions (1)

``` r
summary(model.lm.3 <- lm(stfdem ~ factor(sex) * lrscale, data = dt))
```
This is the same as the more verbose syntax using `:`.

``` r
summary(model.lm.4 <- lm(stfdem ~ factor(sex) + lrscale + factor(sex):lrscale, 
                         data = dt))
```

But to get *full* marginal effects, use the `f1 / f2` interaction syntax:

``` r
summary(model.lm.5 <- lm(stfdem ~ factor(sex) / lrscale, data = dt))
```

---
# Models with interactions (2)
<svg viewBox="0 0 576 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M550.5 241l-50.089-86.786c1.071-2.142 1.875-4.553 1.875-7.232 0-8.036-6.696-14.733-14.732-15.001l-55.447-95.893c.536-1.607 1.071-3.214 1.071-4.821 0-8.571-6.964-15.268-15.268-15.268-4.821 0-8.839 2.143-11.786 5.625H299.518C296.839 18.143 292.821 16 288 16s-8.839 2.143-11.518 5.625H170.411C167.464 18.143 163.447 16 158.625 16c-8.303 0-15.268 6.696-15.268 15.268 0 1.607.536 3.482 1.072 4.821l-55.983 97.233c-5.356 2.41-9.107 7.5-9.107 13.661 0 .535.268 1.071.268 1.607l-53.304 92.143c-7.232 1.339-12.59 7.5-12.59 15 0 7.232 5.089 13.393 12.054 15l55.179 95.358c-.536 1.607-.804 2.946-.804 4.821 0 7.232 5.089 13.393 12.054 14.732l51.697 89.732c-.536 1.607-1.071 3.482-1.071 5.357 0 8.571 6.964 15.268 15.268 15.268 4.821 0 8.839-2.143 11.518-5.357h106.875C279.161 493.857 283.447 496 288 496s8.839-2.143 11.518-5.357h107.143c2.678 2.946 6.696 4.821 10.982 4.821 8.571 0 15.268-6.964 15.268-15.268 0-1.607-.267-2.946-.803-4.285l51.697-90.268c6.964-1.339 12.054-7.5 12.054-14.732 0-1.607-.268-3.214-.804-4.821l54.911-95.358c6.964-1.339 12.322-7.5 12.322-15-.002-7.232-5.092-13.393-11.788-14.732zM153.535 450.732l-43.66-75.803h43.66v75.803zm0-83.839h-43.66c-.268-1.071-.804-2.142-1.339-3.214l44.999-47.41v50.624zm0-62.411l-50.357 53.304c-1.339-.536-2.679-1.34-4.018-1.607L43.447 259.75c.535-1.339.535-2.679.535-4.018s0-2.41-.268-3.482l51.965-90c2.679-.268 5.357-1.072 7.768-2.679l50.089 51.965v92.946zm0-102.322l-45.803-47.41c1.339-2.143 2.143-4.821 2.143-7.767 0-.268-.268-.804-.268-1.072l43.928-15.804v72.053zm0-80.625l-43.66 15.804 43.66-75.536v59.732zm326.519 39.108l.804 1.339L445.5 329.125l-63.75-67.232 98.036-101.518.268.268zM291.75 355.107l11.518 11.786H280.5l11.25-11.786zm-.268-11.25l-83.303-85.446 79.553-84.375 83.036 87.589-79.286 82.232zm5.357 5.893l79.286-82.232 67.5 71.25-5.892 28.125H313.714l-16.875-17.143zM410.411 44.393c1.071.536 2.142 1.072 3.482 1.34l57.857 100.714v.536c0 2.946.803 5.624 2.143 7.767L376.393 256l-83.035-87.589L410.411 44.393zm-9.107-2.143L287.732 162.518l-57.054-60.268 166.339-60h4.287zm-123.483 0c2.678 2.678 6.16 4.285 10.179 4.285s7.5-1.607 10.179-4.285h75L224.786 95.821 173.893 42.25h103.928zm-116.249 5.625l1.071-2.142a33.834 33.834 0 0 0 2.679-.804l51.161 53.84-54.911 19.821V47.875zm0 79.286l60.803-21.964 59.732 63.214-79.553 84.107-40.982-42.053v-83.304zm0 92.678L198 257.607l-36.428 38.304v-76.072zm0 87.858l42.053-44.464 82.768 85.982-17.143 17.678H161.572v-59.196zm6.964 162.053c-1.607-1.607-3.482-2.678-5.893-3.482l-1.071-1.607v-89.732h99.91l-91.607 94.821h-1.339zm129.911 0c-2.679-2.41-6.428-4.285-10.447-4.285s-7.767 1.875-10.447 4.285h-96.429l91.607-94.821h38.304l91.607 94.821H298.447zm120-11.786l-4.286 7.5c-1.339.268-2.41.803-3.482 1.339l-89.196-91.875h114.376l-17.412 83.036zm12.856-22.232l12.858-60.803h21.964l-34.822 60.803zm34.822-68.839h-20.357l4.553-21.16 17.143 18.214c-.535.803-1.071 1.874-1.339 2.946zm66.161-107.411l-55.447 96.697c-1.339.535-2.679 1.071-4.018 1.874l-20.625-21.964 34.554-163.928 45.803 79.286c-.267 1.339-.803 2.678-.803 4.285 0 1.339.268 2.411.536 3.75z"></path></svg> The models are equivalent, but in the third case we do not have to compute ourselves the effect of *Left/Right* for different values of *Gender*, and we get the correct standard errors and p-values computed for us directly for every stratum of *Gender*. This also makes it easier to plot meaningful marginal effects.

In terms of syntax, `f1 / x2` is equivalent to `f1 + f1:x2`, while `f1 * x2` is equivalent to `f1 + x2 + f1:x2`. So `f1 / x2` drops the parent effect of `x2` from the formula.

For details and extensions, see [this post](https://grantmcdermott.com/posts/interaction-effects/) by Grant McDermott. The general issue has to do with the difference between crossed and nested interactions (see  [here](https://stackoverflow.com/questions/32616762/defining-an-infix-operator-for-use-within-a-formula/32682826#32682826)).

---
# Tables of model results (1)
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M464 32H48C21.49 32 0 53.49 0 80v352c0 26.51 21.49 48 48 48h416c26.51 0 48-21.49 48-48V80c0-26.51-21.49-48-48-48zM224 416H64v-96h160v96zm0-160H64v-96h160v96zm224 160H288v-96h160v96zm0-160H288v-96h160v96z"></path></svg> The `summary()` function prints tables of results in the console. To export model results outside `R`, we can use several different packages. It is straightforward to export directly into `latex`, `markdown` and `html` formats. There is also a way to get the regression tables output into everyone's favorite word processor, **MSWord**. With [`modelsummary`](https://modelsummary.com/vignettes/modelsummary.html) or with [`gtsummary`](https://www.danieldsjoberg.com/gtsummary/index.html), we can directly make a `docx` file. There is an alternative that is just a little bit more circumvent: 
- first, we (install and) load the [`stargazer`]((https://cran.r-project.org/web/packages/stargazer/vignettes/stargazer.pdf)) package (or `modelsummary` or `gtsummary`)
- second, we generate a table of results as a small `html` file 
- third, we locate the html file and open it in any browser
- fourth, we copy and paste the html table into MSWord 
- fifth, we are done after some optional tinkering with the formatting

---
# Tables of model results (2)
If convenience wasn't enough, these packages come with other perks, such as computing robust standard errors on the fly, or exponentiating coefficients from logistic regression models. Let's summarize the results of our model of voting behavior:

``` r
library(gtsummary)
tbl_regression(model.glm.1, exponentiate = TRUE)
```

<div id="mgrctzpjta" style="padding-left:0px;padding-right:0px;padding-top:10px;padding-bottom:10px;overflow-x:auto;overflow-y:auto;width:auto;height:auto;">
<style>#mgrctzpjta table {
  font-family: system-ui, 'Segoe UI', Roboto, Helvetica, Arial, sans-serif, 'Apple Color Emoji', 'Segoe UI Emoji', 'Segoe UI Symbol', 'Noto Color Emoji';
  -webkit-font-smoothing: antialiased;
  -moz-osx-font-smoothing: grayscale;
}

#mgrctzpjta thead, #mgrctzpjta tbody, #mgrctzpjta tfoot, #mgrctzpjta tr, #mgrctzpjta td, #mgrctzpjta th {
  border-style: none;
}

#mgrctzpjta p {
  margin: 0;
  padding: 0;
}

#mgrctzpjta .gt_table {
  display: table;
  border-collapse: collapse;
  line-height: normal;
  margin-left: auto;
  margin-right: auto;
  color: #333333;
  font-size: 16px;
  font-weight: normal;
  font-style: normal;
  background-color: #FFFFFF;
  width: auto;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #A8A8A8;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #A8A8A8;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
}

#mgrctzpjta .gt_caption {
  padding-top: 4px;
  padding-bottom: 4px;
}

#mgrctzpjta .gt_title {
  color: #333333;
  font-size: 125%;
  font-weight: initial;
  padding-top: 4px;
  padding-bottom: 4px;
  padding-left: 5px;
  padding-right: 5px;
  border-bottom-color: #FFFFFF;
  border-bottom-width: 0;
}

#mgrctzpjta .gt_subtitle {
  color: #333333;
  font-size: 85%;
  font-weight: initial;
  padding-top: 3px;
  padding-bottom: 5px;
  padding-left: 5px;
  padding-right: 5px;
  border-top-color: #FFFFFF;
  border-top-width: 0;
}

#mgrctzpjta .gt_heading {
  background-color: #FFFFFF;
  text-align: center;
  border-bottom-color: #FFFFFF;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
}

#mgrctzpjta .gt_bottom_border {
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#mgrctzpjta .gt_col_headings {
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
}

#mgrctzpjta .gt_col_heading {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: normal;
  text-transform: inherit;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: bottom;
  padding-top: 5px;
  padding-bottom: 6px;
  padding-left: 5px;
  padding-right: 5px;
  overflow-x: hidden;
}

#mgrctzpjta .gt_column_spanner_outer {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: normal;
  text-transform: inherit;
  padding-top: 0;
  padding-bottom: 0;
  padding-left: 4px;
  padding-right: 4px;
}

#mgrctzpjta .gt_column_spanner_outer:first-child {
  padding-left: 0;
}

#mgrctzpjta .gt_column_spanner_outer:last-child {
  padding-right: 0;
}

#mgrctzpjta .gt_column_spanner {
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  vertical-align: bottom;
  padding-top: 5px;
  padding-bottom: 5px;
  overflow-x: hidden;
  display: inline-block;
  width: 100%;
}

#mgrctzpjta .gt_spanner_row {
  border-bottom-style: hidden;
}

#mgrctzpjta .gt_group_heading {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: middle;
  text-align: left;
}

#mgrctzpjta .gt_empty_group_heading {
  padding: 0.5px;
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  vertical-align: middle;
}

#mgrctzpjta .gt_from_md > :first-child {
  margin-top: 0;
}

#mgrctzpjta .gt_from_md > :last-child {
  margin-bottom: 0;
}

#mgrctzpjta .gt_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  margin: 10px;
  border-top-style: solid;
  border-top-width: 1px;
  border-top-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: middle;
  overflow-x: hidden;
}

#mgrctzpjta .gt_stub {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-right-style: solid;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  padding-left: 5px;
  padding-right: 5px;
}

#mgrctzpjta .gt_stub_row_group {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-right-style: solid;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  padding-left: 5px;
  padding-right: 5px;
  vertical-align: top;
}

#mgrctzpjta .gt_row_group_first td {
  border-top-width: 2px;
}

#mgrctzpjta .gt_row_group_first th {
  border-top-width: 2px;
}

#mgrctzpjta .gt_summary_row {
  color: #333333;
  background-color: #FFFFFF;
  text-transform: inherit;
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
}

#mgrctzpjta .gt_first_summary_row {
  border-top-style: solid;
  border-top-color: #D3D3D3;
}

#mgrctzpjta .gt_first_summary_row.thick {
  border-top-width: 2px;
}

#mgrctzpjta .gt_last_summary_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#mgrctzpjta .gt_grand_summary_row {
  color: #333333;
  background-color: #FFFFFF;
  text-transform: inherit;
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
}

#mgrctzpjta .gt_first_grand_summary_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-top-style: double;
  border-top-width: 6px;
  border-top-color: #D3D3D3;
}

#mgrctzpjta .gt_last_grand_summary_row_top {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-bottom-style: double;
  border-bottom-width: 6px;
  border-bottom-color: #D3D3D3;
}

#mgrctzpjta .gt_striped {
  background-color: rgba(128, 128, 128, 0.05);
}

#mgrctzpjta .gt_table_body {
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#mgrctzpjta .gt_footnotes {
  color: #333333;
  background-color: #FFFFFF;
  border-bottom-style: none;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
}

#mgrctzpjta .gt_footnote {
  margin: 0px;
  font-size: 90%;
  padding-top: 4px;
  padding-bottom: 4px;
  padding-left: 5px;
  padding-right: 5px;
}

#mgrctzpjta .gt_sourcenotes {
  color: #333333;
  background-color: #FFFFFF;
  border-bottom-style: none;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
}

#mgrctzpjta .gt_sourcenote {
  font-size: 90%;
  padding-top: 4px;
  padding-bottom: 4px;
  padding-left: 5px;
  padding-right: 5px;
}

#mgrctzpjta .gt_left {
  text-align: left;
}

#mgrctzpjta .gt_center {
  text-align: center;
}

#mgrctzpjta .gt_right {
  text-align: right;
  font-variant-numeric: tabular-nums;
}

#mgrctzpjta .gt_font_normal {
  font-weight: normal;
}

#mgrctzpjta .gt_font_bold {
  font-weight: bold;
}

#mgrctzpjta .gt_font_italic {
  font-style: italic;
}

#mgrctzpjta .gt_super {
  font-size: 65%;
}

#mgrctzpjta .gt_footnote_marks {
  font-size: 75%;
  vertical-align: 0.4em;
  position: initial;
}

#mgrctzpjta .gt_asterisk {
  font-size: 100%;
  vertical-align: 0;
}

#mgrctzpjta .gt_indent_1 {
  text-indent: 5px;
}

#mgrctzpjta .gt_indent_2 {
  text-indent: 10px;
}

#mgrctzpjta .gt_indent_3 {
  text-indent: 15px;
}

#mgrctzpjta .gt_indent_4 {
  text-indent: 20px;
}

#mgrctzpjta .gt_indent_5 {
  text-indent: 25px;
}

#mgrctzpjta .katex-display {
  display: inline-flex !important;
  margin-bottom: 0.75em !important;
}

#mgrctzpjta div.Reactable > div.rt-table > div.rt-thead > div.rt-tr.rt-tr-group-header > div.rt-th-group:after {
  height: 0px !important;
}
</style>
<table class="gt_table" data-quarto-disable-processing="false" data-quarto-bootstrap="false">
  <thead>
    <tr class="gt_col_headings">
      <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1" scope="col" id="label"><span class='gt_from_md'><strong>Characteristic</strong></span></th>
      <th class="gt_col_heading gt_columns_bottom_border gt_center" rowspan="1" colspan="1" scope="col" id="estimate"><span class='gt_from_md'><strong>OR</strong></span></th>
      <th class="gt_col_heading gt_columns_bottom_border gt_center" rowspan="1" colspan="1" scope="col" id="conf.low"><span class='gt_from_md'><strong>95% CI</strong></span></th>
      <th class="gt_col_heading gt_columns_bottom_border gt_center" rowspan="1" colspan="1" scope="col" id="p.value"><span class='gt_from_md'><strong>p-value</strong></span></th>
    </tr>
  </thead>
  <tbody class="gt_table_body">
    <tr><td headers="label" class="gt_row gt_left">How satisfied with the way democracy works in country</td>
<td headers="estimate" class="gt_row gt_center">0.98</td>
<td headers="conf.low" class="gt_row gt_center">0.97, 1.00</td>
<td headers="p.value" class="gt_row gt_center">0.017</td></tr>
    <tr><td headers="label" class="gt_row gt_left">sex</td>
<td headers="estimate" class="gt_row gt_center"><br /></td>
<td headers="conf.low" class="gt_row gt_center"><br /></td>
<td headers="p.value" class="gt_row gt_center"><br /></td></tr>
    <tr><td headers="label" class="gt_row gt_left">    female</td>
<td headers="estimate" class="gt_row gt_center">—</td>
<td headers="conf.low" class="gt_row gt_center">—</td>
<td headers="p.value" class="gt_row gt_center"><br /></td></tr>
    <tr><td headers="label" class="gt_row gt_left">    male</td>
<td headers="estimate" class="gt_row gt_center">1.02</td>
<td headers="conf.low" class="gt_row gt_center">0.97, 1.07</td>
<td headers="p.value" class="gt_row gt_center">0.5</td></tr>
    <tr><td headers="label" class="gt_row gt_left">Placement on left right scale</td>
<td headers="estimate" class="gt_row gt_center">1.02</td>
<td headers="conf.low" class="gt_row gt_center">1.01, 1.03</td>
<td headers="p.value" class="gt_row gt_center"><0.001</td></tr>
    <tr><td headers="label" class="gt_row gt_left">Trust in country's parliament</td>
<td headers="estimate" class="gt_row gt_center">1.07</td>
<td headers="conf.low" class="gt_row gt_center">1.06, 1.08</td>
<td headers="p.value" class="gt_row gt_center"><0.001</td></tr>
  </tbody>
  <tfoot class="gt_sourcenotes">
    <tr>
      <td class="gt_sourcenote" colspan="4"><span class='gt_from_md'>Abbreviations: CI = Confidence Interval, OR = Odds Ratio</span></td>
    </tr>
  </tfoot>
  
</table>
</div>

---
# Let's see how to export regression tables with `stargazer`

``` r
m1a<-lm(stfdem ~ lrscale + trstprl + etfruit, data=dt)
m1b<-lm(stfdem ~ lrscale + trstprl + etfruit + sex, data=dt)
m1c<-update (m1b, .~. + stflife)

library(stargazer)
stargazer(m1a, m1b, m1c, type="html", align=T, omit.stat=c("LL","ser","f"), 
          no.space=T, title = 'Models of something very important',
          dep.var.labels=c("Satisfaction wtih democracy"), 
          covariate.labels=c("Left-Right","Trust (nat.)", "Eats fruit", 'Gender (male)', 'Life satisfaction'),
          column.labels=c('Model 1.1','Model 1.2','Model 1.3'), 
          digits=2, initial.zero =T, model.numbers=F, 
          out="./models/model_set1_sg.html")
## 
## <table style="text-align:center"><caption><strong>Models of something very important</strong></caption>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="3"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="3" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="3">Satisfaction wtih democracy</td></tr>
## <tr><td style="text-align:left"></td><td>Model 1.1</td><td>Model 1.2</td><td>Model 1.3</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Left-Right</td><td>0.09<sup>***</sup></td><td>0.09<sup>***</sup></td><td>0.07<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.005)</td><td>(0.005)</td><td>(0.005)</td></tr>
## <tr><td style="text-align:left">Trust (nat.)</td><td>0.59<sup>***</sup></td><td>0.59<sup>***</sup></td><td>0.55<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.004)</td><td>(0.004)</td><td>(0.004)</td></tr>
## <tr><td style="text-align:left">Eats fruit</td><td>-0.04<sup>***</sup></td><td>-0.04<sup>***</sup></td><td>-0.02<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td><td>(0.01)</td><td>(0.01)</td></tr>
## <tr><td style="text-align:left">Gender (male)</td><td></td><td>0.07<sup>***</sup></td><td>0.06<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.02)</td><td>(0.02)</td></tr>
## <tr><td style="text-align:left">Life satisfaction</td><td></td><td></td><td>0.20<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left">Constant</td><td>2.20<sup>***</sup></td><td>2.19<sup>***</sup></td><td>0.91<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.04)</td><td>(0.04)</td><td>(0.05)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>34,418</td><td>34,418</td><td>34,361</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.38</td><td>0.38</td><td>0.40</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.38</td><td>0.38</td><td>0.40</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="3" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
```

---
# Now let's take `modelsummary` for a spin

``` r
library(modelsummary)
models1 = dvnames(list (m1a, m1b, m1c))
modelsummary(models1,  fmt = 2, 
             estimate = "{estimate} [{conf.low}, {conf.high}] p={p.value} {stars}",
             statistic = NULL, 
             title = '',
             #coef_omit = ("cntry"),
             gof_omit = 'DF|Deviance|AIC|BIC|RMSE',
             output="./models/model_set1_ms.docx")
```

---
# Mediation analysis

We can use state-of-the-art causal mediation analysis with the `mediation` package: [details](https://cran.r-project.org/web/packages/mediation/vignettes/mediation.pdf).
- first we specify a model of the response
- then we specify a model of the mediator
- then we let `mediate()` do its magic
- and then we run sensitivity analyses

For an interaction between treatment and mediator just add an interaction in the outcome model: to test the significance of the interaction, use `test.TMint(med.out, conf.level = .95)`. Use `medsens()` for sensitivity analysis.

---
# Mediation analysis in practice (1)
Let's build a model of *Life satisfaction* as a function of *health status* as mediated by *frequency of doing sports*. Scroll down to see the results.

``` r
library(mediation)
dt$health <- 6-dt$health # invert the coding to correspond with the label
dt.s <- dt %>% dplyr::select(health, agea, sex, dosprt, stflife) %>% na.omit()
med.model <- lm(dosprt ~ health + agea + sex, data=dt.s) # model of mediator
out.model <- lm(stflife ~ health + agea + sex + dosprt, data=dt.s) # model of outcome

final.model <- mediate(med.model, out.model, treat = "health", 
                       mediator = "dosprt", robustSE = TRUE, sims = 100) 
summary(final.model)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME             0.0404       0.0362         0.04  <2e-16 ***
## ADE              0.6579       0.6371         0.68  <2e-16 ***
## Total Effect     0.6983       0.6779         0.72  <2e-16 ***
## Prop. Mediated   0.0576       0.0518         0.06  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 39164 
## 
## 
## Simulations: 100
```
---
# Mediation analysis in practice (2)
Plotting the results is easy:

.pull-left[

``` r
plot(final.model, xlim=c(0,1),
     main='Effect of Health on 
     Life Satisfaction\nDirect and 
     Mediated by Doing Sports')
```
]

.pull-right[
<img src="data:image/png;base64,#R-Tutorial-Session-3_files/figure-html/unnamed-chunk-26-1.png" width="504" />
]
---
# More modelling?
.center[![trachoma](data:image/png;base64,#figs/zoolander1.jpg)]

---
# Multilevel models
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M12.41 148.02l232.94 105.67c6.8 3.09 14.49 3.09 21.29 0l232.94-105.67c16.55-7.51 16.55-32.52 0-40.03L266.65 2.31a25.607 25.607 0 0 0-21.29 0L12.41 107.98c-16.55 7.51-16.55 32.53 0 40.04zm487.18 88.28l-58.09-26.33-161.64 73.27c-7.56 3.43-15.59 5.17-23.86 5.17s-16.29-1.74-23.86-5.17L70.51 209.97l-58.1 26.33c-16.55 7.5-16.55 32.5 0 40l232.94 105.59c6.8 3.08 14.49 3.08 21.29 0L499.59 276.3c16.55-7.5 16.55-32.5 0-40zm0 127.8l-57.87-26.23-161.86 73.37c-7.56 3.43-15.59 5.17-23.86 5.17s-16.29-1.74-23.86-5.17L70.29 337.87 12.41 364.1c-16.55 7.5-16.55 32.5 0 40l232.94 105.59c6.8 3.08 14.49 3.08 21.29 0L499.59 404.1c16.55-7.5 16.55-32.5 0-40z"></path></svg> We can let our modeling span multiple levels with the `lme4` package.

``` r
library(lme4)
library(arm)

# varying intercept with individual level predictor
ml0 <- lmer(stfdem ~ lrscale + (1 | cntry), data=dt)

# varying intercept, varying slope
ml1 <- lmer(stfdem ~ lrscale + (1 + lrscale | cntry), data=dt) 
```

---
# America's next top (multilevel) model
We can access the estimated coefficients, as well as the computed fixed effects and random effects. For mode details, see the excellent book by [Gelman and Hill (2007)](https://amzn.to/2Hgz5fA). For a great intuitive introduction to multilevel models, see this [site](http://mfviz.com/hierarchical-models/) by Michael Freeman.

``` r
summary(ml1)
display(ml1)
coef(ml1) # estimated coefficients
fixef(ml1) # fixed effects
ranef(ml1) # random effects
```

---
# Understanding multilevel models
To understand, illustrate and export results from multilevel models, use the `sjPlot` package and the `tab_model()` function in particular. For details, see [here](https://strengejacke.github.io/sjPlot/articles/tab_mixed.html).

``` r
library(sjPlot)
tab_model(ml1)
```

<table style="border-collapse:collapse; border:none;">
<tr>
<th style="border-top: double; text-align:center; font-style:normal; font-weight:bold; padding:0.2cm;  text-align:left; ">&nbsp;</th>
<th colspan="3" style="border-top: double; text-align:center; font-style:normal; font-weight:bold; padding:0.2cm; ">How satisfied with the<br>way democracy works in<br>country</th>
</tr>
<tr>
<td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal;  text-align:left; ">Predictors</td>
<td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal;  ">Estimates</td>
<td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal;  ">CI</td>
<td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal;  ">p</td>
</tr>
<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">(Intercept)</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center;  ">4.78</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center;  ">4.10&nbsp;&ndash;&nbsp;5.45</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center;  "><strong>&lt;0.001</strong></td>
</tr>
<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">Placement on left right<br>scale</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center;  ">0.10</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center;  ">0.02&nbsp;&ndash;&nbsp;0.18</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center;  "><strong>0.014</strong></td>
</tr>
<tr>
<td colspan="4" style="font-weight:bold; text-align:left; padding-top:.8em;">Random Effects</td>
</tr>

<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">&sigma;<sup>2</sup></td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">5.19</td>
</tr>

<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">&tau;<sub>00</sub> <sub>cntry</sub></td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">2.84</td>

<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">&tau;<sub>11</sub> <sub>cntry.lrscale</sub></td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">0.04</td>

<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">&rho;<sub>01</sub> <sub>cntry</sub></td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">-0.85</td>

<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">N <sub>cntry</sub></td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">24</td>
<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm; border-top:1px solid;">Observations</td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left; border-top:1px solid;" colspan="3">34733</td>
</tr>
<tr>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; padding-top:0.1cm; padding-bottom:0.1cm;">Marginal R<sup>2</sup> / Conditional R<sup>2</sup></td>
<td style=" padding:0.2cm; text-align:left; vertical-align:top; padding-top:0.1cm; padding-bottom:0.1cm; text-align:left;" colspan="3">0.008 / 0.191</td>
</tr>

</table>
---
# Analysis of randomized experiments (1)
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M504.971 359.029c9.373 9.373 9.373 24.569 0 33.941l-80 79.984c-15.01 15.01-40.971 4.49-40.971-16.971V416h-58.785a12.004 12.004 0 0 1-8.773-3.812l-70.556-75.596 53.333-57.143L352 336h32v-39.981c0-21.438 25.943-31.998 40.971-16.971l80 79.981zM12 176h84l52.781 56.551 53.333-57.143-70.556-75.596A11.999 11.999 0 0 0 122.785 96H12c-6.627 0-12 5.373-12 12v56c0 6.627 5.373 12 12 12zm372 0v39.984c0 21.46 25.961 31.98 40.971 16.971l80-79.984c9.373-9.373 9.373-24.569 0-33.941l-80-79.981C409.943 24.021 384 34.582 384 56.019V96h-58.785a12.004 12.004 0 0 0-8.773 3.812L96 336H12c-6.627 0-12 5.373-12 12v56c0 6.627 5.373 12 12 12h110.785c3.326 0 6.503-1.381 8.773-3.812L352 176h32z"></path></svg> As you would expect, the analysis of experiments is very easy. You can of course use `lm()` to retrieve average treatment effects. If you want to run analysis of variance, the function `stats::anova()` in `base R` takes a fitted model object (so it takes the output of `lm()`, for example) with **Type I** sums of squares.

---
# Analysis of randomized experiments (2)
It gets a tad more complicated when we want to include covariates. The [recommended approach](https://alexandercoppock.com/Green-Lab-SOP/Green_Lab_SOP.html) is to specify a model wtih robust standard errors in which the treatment is interacted with the demeaned covariates. This leads to the so-called [*Lin regression*](http://web.archive.org/web/20150505184132/http://blogs.worldbank.org/impactevaluations/node/847.). You can implement this approach via the `estimatr` package and its function `lm_lin()`. Alternatively, use `lm_robust()`, but don't forget to demean the covariates.

``` r
library(estimatr)
# Hypothetical experimental dataset `d`
summary(em1a<-lm_lin(dv.index ~ condition, covariates = ~ age, data=d)) # apprach I

d$age_demeaned = scale(d$age, scale=FALSE) # to demean (center) the covariate
summary(em1b<-lm_robust(dv.index ~ condition * age_demeaned, data=d)) # approach II
```

---
# Analysis of randomized experiments (3)
Before analysis comes design. The [`DeclareDesign`](https://declaredesign.org/r/declaredesign/) package supports the design of complex experiments, including statistical power analyses, decisions about blocked randomization, spillover considerations, etc. It is a very powerful tool, but it is not very easy to master. Try it out the next time you plan your experiments.

Conjoint experimental designs are all the rage these days. Personally, I am not a fan, but if you want to learn how to analyse data from conjoint experiments in R, follow [this tutorial in two parts](https://www.andrewheiss.com/blog/2023/08/12/conjoint-multilevel-multinomial-guide/) by Andrew Heiss.

---
# Structural Equations Models (SEM)
<svg viewBox="0 0 640 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M384 320H256c-17.67 0-32 14.33-32 32v128c0 17.67 14.33 32 32 32h128c17.67 0 32-14.33 32-32V352c0-17.67-14.33-32-32-32zM192 32c0-17.67-14.33-32-32-32H32C14.33 0 0 14.33 0 32v128c0 17.67 14.33 32 32 32h95.72l73.16 128.04C211.98 300.98 232.4 288 256 288h.28L192 175.51V128h224V64H192V32zM608 0H480c-17.67 0-32 14.33-32 32v128c0 17.67 14.33 32 32 32h128c17.67 0 32-14.33 32-32V32c0-17.67-14.33-32-32-32z"></path></svg> There are (at least) two packages for doing SEM: `sem` and `lavaan`. More on the SEM package [here](https://socialsciences.mcmaster.ca/jfox/Misc/sem/SEM-paper.pdf).
I recommend working with [`lavaan`](https://www.jstatsoft.org/article/view/v048i02).
.pull-left[

``` r
library(lavaan)
l.model <- '
   # latent variables
     lat1 =~ trstun + trstep + trstprl
     lat2 =~ stfdem + lrscale
   # regressions
     lat2 ~ lat1
   # residual covariances
     lrscale ~~ trstun
'
```
]
.pull-right[

``` r
l.fit <- lavaan::sem(l.model, data=dt)
#summary(l.fit)
parameterEstimates(l.fit)
##        lhs op     rhs        est se  z pvalue ci.lower ci.upper
## 1     lat1 =~  trstun      1.000  0 NA     NA        1        1
## 2     lat1 =~  trstep      1.046 NA NA     NA       NA       NA
## 3     lat1 =~ trstprl      0.836 NA NA     NA       NA       NA
## 4     lat2 =~  stfdem      1.000  0 NA     NA        1        1
## 5     lat2 =~ lrscale      0.000 NA NA     NA       NA       NA
## 6     lat2  ~    lat1      0.688 NA NA     NA       NA       NA
## 7   trstun ~~ lrscale      0.082 NA NA     NA       NA       NA
## 8   trstun ~~  trstun      2.103 NA NA     NA       NA       NA
## 9   trstep ~~  trstep      1.544 NA NA     NA       NA       NA
## 10 trstprl ~~ trstprl      3.741 NA NA     NA       NA       NA
## 11  stfdem ~~  stfdem -23643.118 NA NA     NA       NA       NA
## 12 lrscale ~~ lrscale      5.038 NA NA     NA       NA       NA
## 13    lat1 ~~    lat1      4.317 NA NA     NA       NA       NA
## 14    lat2 ~~    lat2  23647.437 NA NA     NA       NA       NA
```
]

---
# PCA and factor analysis
R has extensive capabilities for factor analysis, PCA and other dimension-reduction techniques, such as item response models.

``` r
library(psych)
# principal components analysis
pca(dt.com)

# factor analysis
model.fa1<-fa(dt.com, nfactors=2, rotate='varimax', scores='regression') 
print(model.fa1)
model.fa1$scores[1:5,]
```
<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M496 448H16c-8.84 0-16 7.16-16 16v32c0 8.84 7.16 16 16 16h480c8.84 0 16-7.16 16-16v-32c0-8.84-7.16-16-16-16zm-304-64l-64-32 64-32 32-64 32 64 64 32-64 32-16 32h208l-86.41-201.63a63.955 63.955 0 0 1-1.89-45.45L416 0 228.42 107.19a127.989 127.989 0 0 0-53.46 59.15L64 416h144l-16-32zm64-224l16-32 16 32 32 16-32 16-16 32-16-32-32-16 32-16z"></path></svg> **Protip:** Use tetrachronic correlations for variables with few levels.

---
# Yet more models
If what we reviewed above doesn't have you covered, consider:
- models for categorical data (ordered and unordered) with the `MASS` package
- event history models (survival analysis) with the `survival` package
- instrumental variables and other econometric models with the `AER` package
- generalized additive models with `gam`
- survey design with the `survey` package

---
# Don't forget to save your work
<svg viewBox="0 0 448 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M433.941 129.941l-83.882-83.882A48 48 0 0 0 316.118 32H48C21.49 32 0 53.49 0 80v352c0 26.51 21.49 48 48 48h352c26.51 0 48-21.49 48-48V163.882a48 48 0 0 0-14.059-33.941zM272 80v80H144V80h128zm122 352H54a6 6 0 0 1-6-6V86a6 6 0 0 1 6-6h42v104c0 13.255 10.745 24 24 24h176c13.255 0 24-10.745 24-24V83.882l78.243 78.243a6 6 0 0 1 1.757 4.243V426a6 6 0 0 1-6 6zM224 232c-48.523 0-88 39.477-88 88s39.477 88 88 88 88-39.477 88-88-39.477-88-88-88zm0 128c-22.056 0-40-17.944-40-40s17.944-40 40-40 40 17.944 40 40-17.944 40-40 40z"></path></svg> Save code you've written and wanna use in the future.
<svg viewBox="0 0 448 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M433.941 129.941l-83.882-83.882A48 48 0 0 0 316.118 32H48C21.49 32 0 53.49 0 80v352c0 26.51 21.49 48 48 48h352c26.51 0 48-21.49 48-48V163.882a48 48 0 0 0-14.059-33.941zM272 80v80H144V80h128zm122 352H54a6 6 0 0 1-6-6V86a6 6 0 0 1 6-6h42v104c0 13.255 10.745 24 24 24h176c13.255 0 24-10.745 24-24V83.882l78.243 78.243a6 6 0 0 1 1.757 4.243V426a6 6 0 0 1-6 6zM224 232c-48.523 0-88 39.477-88 88s39.477 88 88 88 88-39.477 88-88-39.477-88-88-88zm0 128c-22.056 0-40-17.944-40-40s17.944-40 40-40 40 17.944 40 40-17.944 40-40 40z"></path></svg> Save data that you've created or modified.
<svg viewBox="0 0 448 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg">  <path d="M433.941 129.941l-83.882-83.882A48 48 0 0 0 316.118 32H48C21.49 32 0 53.49 0 80v352c0 26.51 21.49 48 48 48h352c26.51 0 48-21.49 48-48V163.882a48 48 0 0 0-14.059-33.941zM272 80v80H144V80h128zm122 352H54a6 6 0 0 1-6-6V86a6 6 0 0 1 6-6h42v104c0 13.255 10.745 24 24 24h176c13.255 0 24-10.745 24-24V83.882l78.243 78.243a6 6 0 0 1 1.757 4.243V426a6 6 0 0 1-6 6zM224 232c-48.523 0-88 39.477-88 88s39.477 88 88 88 88-39.477 88-88-39.477-88-88-88zm0 128c-22.056 0-40-17.944-40-40s17.944-40 40-40 40 17.944 40 40-17.944 40-40 40z"></path></svg> Save graphs, tables and model output.

---
# How to get in touch?