Quality of Government dataset

# the standard, cross-section dataset
# qog_path <- "http://www.qogdata.pol.gu.se/data/qog_std_cs_jan19.dta"
qog <- rio::import(qog_path)

Summarise the data

qog_summary <- sjmisc::descr(qog)

1. Simulate the data many, many times.

2. What percentage of these simulations lead to the same regression line as the one from the real data?

Answer to question 2 is the p-value.

Common misinterpretation:

• The p-value is the probability that the null hypothesis is false.
• A small enough p-value means that we can be certain that the relationship between the variables is non-zero.

\[ p(h|I, D) = p(h|I) \frac{p(D|h,I)}{p(D|I)} \]

• \(h\): the hypothesis
• \(D\): the empirical data
• \(I\): the background information we had prior to analysing the data

\[ p(h|I, D) = p(h|I) \frac{p(D|h,I)}{p(D|I)} \]

• Posterior probability: \(p(h|I, D)\)
• Prior probability: \(p(h|I)\)
• Likelihood: \(p(D|h,I)\)
• Credibility of the data: \(p(D|I)\)

\[ p(h|I, D) = p(h|I) \frac{\color{red}{p(D|h,I)}}{p(D|I)} \]

• Posterior probability: \(p(h|I, D)\)
• Prior probability: \(p(h|I)\)
• Likelihood: \(p(D|h,I) \longleftarrow\) This is the p-value
• Credibility of the data: \(p(D|I)\)

\[ p(h|I, D) = p(h|I) \frac{p(D|h,I)}{p(D|I)} \] Posterior prob = p-value, only when we don’t have an informative prior

• Why report the p-value intead of the posterior probability?

• Codifying prior knowkedge into a numeric probability is HARD. Calculating the likelihood (the p-value) is relatively easy.

• Consequence: Sometimes we don’t find the regression results credible due to our prior probability (not included in the statistical analysis)

• Most common case: Evaluate the sign of the coefficients – do they make some intuitive sense?

reg1 <- qog %>% lm(unna_gdppc ~ fi_legprop_pd, data = .)
tab_model(reg1, show.se = TRUE)

• v is the omitted variable;
• notice the 1.5 proportionality relation between y and x

v <- rnorm(1000)
x <- v + rnorm(1000)
y <- 1.5 * x + v + rnorm(1000)

Function rnorm generates random numbers from a normal distribution

• r: random
• norm: normal distribution

plot(x, y)

plot(x, v)

plot(v, y)

Notice the coefficient of x is not 1.5

m1 <- lm(y ~ x)
m1

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##    -0.05064      1.99526

Now the coefficient of x is 1.5

m2 <- lm(y ~ x + v)
m2

## 
## Call:
## lm(formula = y ~ x + v)
## 
## Coefficients:
## (Intercept)            x            v  
##    -0.06678      1.51609      0.96514

tab_model(m1, m2)

tab_model(m1, m2, collapse.ci = TRUE)

tab_model(m1, m2, p.style = "asterisk", collapse.ci = TRUE)

tab_model(m1, m2, p.style = "asterisk", 
          show.se = TRUE, show.ci = FALSE, collapse.se = TRUE)

tab_model(m1, m2, p.style = "asterisk", 
          show.se = TRUE, collapse.se = TRUE)

plot_models(m1, m2, show.values = TRUE)

Repeat the analysis a 1000 times to prove that the above omitted variable bias was not an accident.

B <- map_dbl(1:1000, function(i) {
  
    v <- rnorm(1000)
    x <- v + rnorm(1000)
    y <- 1.5 * x + v + rnorm(1000)
    
    lm(y ~ x)$coefficients[[2]]
    
  }) %>% data.frame(beta_x = .)

m1 <- lm(y ~ x)

ggplot(B, aes(x = beta_x)) +
  geom_density() +
  geom_vline(xintercept = mean(B$beta_x))