3.1 Data: Variable Descriptions

• UNSATI - Overall life satisfaction, most recent measure during 2006-2009 (0-10, higher means greater life satisfaction)

• UN_GINI - Distribution of income among individuals or households 2000-2010. 0 represents absolute equality; 100 represents absolute inequality.

• UNPOP10 - Total population in millions, 2010

• UNMAGE10 - Median age in years, 2010

• UNMPI - Multidimensional Poverty Index value, 2000-2008

• UNGDPCAP - Gross domestic product (GDP) per capita in US dollars, 2008

• I_RELIGION - Largest religion by proportion; (1 - Catholic, 2 - Orth. Christian, 3 - Other Christian, 4 - Muslim, 5 - Buddhist, 6 - Other)

• CIAOIL_P - Oil production, in barrels per day

• UNLIT - Adult literacy rate (percentage of population age 15 or older)

• CIALAREA - Land area in square kilometers

• UNPOLFRE - Democracy Score (0 - Authoritarian, 1 - Democratic w/ no alternation, 2 - Democratic)

• UNADJNS - Adjusted net savings

We can quickly look at the geospatial distribution of our dependent variable UNSATI with the help of gvisGeoChart. The graph shows a wide variance in overall life satisfaction distributed across the world, with the most happy people living in Western countries and in some oil producing states like Saudi Arabia. This matches an obvious and strong correlation between poor countries and poor life satisfaction which is evident from most of the African continent.

df_worldplot <- dplyr::select(df, CIANAME, UNSATI) %>% na.omit(df_worldplot) # prepare new df for graph

worldplot <- gvisGeoChart(df_worldplot, locationvar = "CIANAME", # gvis follows certain naming conventions
                          colorvar = "UNSATI", # our dependent variable
                          options=list(width=750, 
                                       height=600))

cat(worldplot$html$chart, file = "worldplot.html") # work around to get our graph to run inside our RMD
htmltools::includeHTML("worldplot.html") # load the graph here

7.1 Analysis: Calculate Slope of Best-Fit Line

We can find the slope apparent in the 2-way scatter above through the following calculation,

\[ \beta = \frac{\sum_{i=1}^N (y_i - \bar{y})(x_i - \bar{x})}{\sum_{i=1}^N (x_i - \bar{x})^2} \]

Whereby:

• \(y_i\) is the life satisfaction of the i-th country in the dataset

• \(\bar{y}\) is the mean life satisfaction

• \(x_i\) is the adult literacy rate of the i-th country in the dataset

• \(\bar{x}\) is the mean adult literacy rate

• \(\sum_{i=1}^N\) is a sum across all observations

Translating this into R, we do the following:

df2 = dplyr::select(df, UNLIT, UNSATI) # select X and Y variables
df2 <- na.omit(df2) # remove NA
df2 <- mutate(df2, mean_literacy = mean(UNLIT),
               mean_satisfaction = mean(UNSATI),
               numerator = (UNLIT - mean_literacy) * (UNSATI - mean_satisfaction),
               denominator = (UNLIT - mean_literacy)^2)
coef <- sum(df2$numerator) / sum(df2$denominator)

Now to check our work, we use the lm function.

# check the calculation above
fit1 <- lm(UNSATI ~ UNLIT, data = df2) # fit a model

And then we compare our by-hand coefficient to that of the lm function.

c(coef, fit1$coefficients[2]) # match

##                 UNLIT 
## 0.03661856 0.03661856

While single linear regression is relatively straight forward, multiple linear regression requires a few more steps and the use of matrix multiplication.

7.2 Analysis: Multivariate Regression by Hand

df_hand <- dplyr::select(df, UNSATI, UNLIT, 
                         UNGDPCAP, UNPOP10, UN_GINI) # pull variables we will use in the model
df_hand <- na.omit(df_hand) # remove NA

Y <- as.matrix(df_hand$UNSATI) # place Y in a vector

X <- data.frame(intercept = 1, 
                       literacy = df_hand$UNLIT, 
                       gdp = df_hand$UNGDPCAP, 
                       pop = df_hand$UNPOP10, 
                       gini = df_hand$UN_GINI) # place X variables in a data frame

X <- data.matrix(X) # convert to matrix

Through linear algebra, the formula for Ordinary Least Squares is parsimoniously listed below,

\[ \hat{\beta} = (X'X)^{-1}X'Y \] Whereby the formula is read as: X prime X inverse \(\times\) X prime Y

• To generate the prime of a matrix, we transpose it. In R, we use the command t(X) to the transpose of the matrix X

• To matrix multiply, we use the command X %*% Y which matrix-multiplies matrix X by matrix Y

• To generate the inverse of a matrix, we use the command solve(X) calculates the matrix-inverse of matrix X

We can calculate multiple linear regression by-hand if we follow the five steps below:

In R, each step looks like the following:

# Step 1 - Transpose X
transpose_X <- t(X)

# Step 2 - Matrix Multiply X' by Y
x_prime_y <- transpose_X %*% Y

# Step 3 - Matrix Multiply X' by X
x_prime_x <- transpose_X %*% X

# Step 4 - X Prime X - Inverse
x_prime_x_inv <- solve(x_prime_x)

# Step 5 - multiply X Prime X Inverse by X Prime Y
betas <- x_prime_x_inv %*% x_prime_y

However, we can use a simplified approach to calculate \(\hat{\beta}\) in one line as shown below:

betas <- solve(t(X) %*% X) %*% t(X) %*% Y
stargazer(betas, type = 'text')

## 
## =================
## intercept  1.517 
## literacy   0.029 
## gdp       0.00005
## pop        0.001 
## gini       0.032 
## -----------------

To see that we did calculation correctly, compare the results to the lm function.

beta_check <- lm(UNSATI ~ UNLIT + UNGDPCAP + UNPOP10 + UN_GINI, data = df)
c(betas, beta_check$coefficients)

##                                                                   (Intercept) 
## 1.516640e+00 2.867659e-02 4.719285e-05 8.055081e-04 3.188312e-02 1.516640e+00 
##        UNLIT     UNGDPCAP      UNPOP10      UN_GINI 
## 2.867659e-02 4.719285e-05 8.055081e-04 3.188312e-02

7.3 Analysis: Interpreting Results

Now that we have a good handle on how we calculate our slope, the next task is its interpretation. Here, we fit a model which we think parsimoniously explains life satisfaction driven by our primary independent variable: literacy. As we know, literacy is a proxy for all sorts of wonderful qualities in adults. Given an ability to read, adults are able to attain better jobs, obtain an education, and importantly, earn more money. While literacy is important, we are also controlling for other factors, such as the natural log of a country’s GDP, the natural log of a country’s population, and the relative distribution of income across households.

fit1 <- lm(UNSATI ~ log(UNLIT) + log(UNGDPCAP) + 
             log(UNPOP10) + UN_GINI, data = df)
summary(fit1)

## 
## Call:
## lm(formula = UNSATI ~ log(UNLIT) + log(UNGDPCAP) + log(UNPOP10) + 
##     UN_GINI, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.41541 -0.53448 -0.00523  0.58033  2.35745 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -4.41417    1.51531  -2.913 0.004442 ** 
## log(UNLIT)     1.32005    0.37985   3.475 0.000765 ***
## log(UNGDPCAP)  0.36334    0.07944   4.573 1.42e-05 ***
## log(UNPOP10)   0.12375    0.07060   1.753 0.082766 .  
## UN_GINI        0.02354    0.01065   2.210 0.029477 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9658 on 97 degrees of freedom
##   (151 observations deleted due to missingness)
## Multiple R-squared:  0.4495, Adjusted R-squared:  0.4268 
## F-statistic:  19.8 on 4 and 97 DF,  p-value: 6.073e-12

stargazer(fit1,
          type = "html",
          title = "Literacy's Relationship with Life Satisfaction",
          dep.var.labels = "Life Satisfaction",
          covariate.labels = c("Literacy", "GDP", "Population", 
                               "Inequality"),
          ci = TRUE)

100*(exp(0.006057) - 1) # convert to percent change in Y

## [1] 0.6075381

We interpret the results, its confidence interval, and the p-value in the following way:

• Effect: A 1% increase in adult literacy is associated with a 0.01% change in life satisfaction, on average, after controlling for GDP, population, and inequality.

• Confidence: If we were to draw repeated random samples and run the same regression on each sample, then 95% of the models’ 95% confidence intervals would contain the true parameter. Therefore, there is a .95 probability that the interval between 0.576 and 2.065 contains the true parameter.

• p-value: Assuming literacy has no effect on life satisfaction, there is a 0.0007 probability that the absolute value of the t statistic would be 3.475 or more.

fit1 <- lm(log(UNSATI) ~ UNLIT + log(UNGDPCAP) + 
             log(UNPOP10) + UN_GINI, data = df)
summary(fit1)

## 
## Call:
## lm(formula = log(UNSATI) ~ UNLIT + log(UNGDPCAP) + log(UNPOP10) + 
##     UN_GINI, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66833 -0.09020  0.00771  0.10900  0.43553 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.534958   0.150942   3.544 0.000608 ***
## UNLIT         0.004710   0.001229   3.832 0.000225 ***
## log(UNGDPCAP) 0.063552   0.016369   3.882 0.000189 ***
## log(UNPOP10)  0.027509   0.013904   1.979 0.050701 .  
## UN_GINI       0.004448   0.002099   2.120 0.036595 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1898 on 97 degrees of freedom
##   (151 observations deleted due to missingness)
## Multiple R-squared:  0.4551, Adjusted R-squared:  0.4327 
## F-statistic: 20.26 on 4 and 97 DF,  p-value: 3.747e-12

100*(exp(0.006057) - 1) # convert to percent change in Y

## [1] 0.6075381

7.4 Analysis: Understanding Control Variables

What do “control” variables actually do in our regressions? The coefficient for our primary independent variable is 0.02308. How did we get this exact number?

First, we will prepare some fresh data for this experiment:

df_controls <- dplyr::select(df, UNSATI, UNLIT, UNGDPCAP, UNPOP10, UN_GINI)
df_controls <- na.omit(df_controls)

Next, we will omit our primary independent variable from the regression.

fit_control1 <- lm(UNSATI ~ log(UNGDPCAP, base = 1.1) + 
             log(UNPOP10) + UN_GINI, data = df_controls)

Then we will save this model’s residuals by:

control_errors1 <- residuals(fit_control1)

Next, we exchange the dependent variable, UNSATI, with our primary IV, UNLIT.

fit_control2 <- lm(UNLIT ~ log(UNGDPCAP, base = 1.1) + 
             log(UNPOP10) + UN_GINI, data = df_controls)

And we save its residuals like so:

control_errors2 <- residuals(fit_control2)

Next, we regress the errors from our first regression on the second. Remember, it is spoken as if “y” is regressed on “x” and not the other way around.

fit_control3 <- lm(control_errors1 ~ control_errors2)

Notice that we now have our beta.

fit_control3$coefficients[2]

## control_errors2 
##      0.02307966

9.1 Analysis: Interactions and Marginal Effects Graphs

plot_ly(data = ME, type = 'scatter', mode = 'lines') %>%
  add_trace(x=~sq_lit, y=~me, name = 'Marginal Effect', line = list(color='rgb(0,100,80)')) %>%
  add_trace(x =~sq_lit, y =~UB, line = list(color = 'transparent'), 
            name = 'Upper CI', showlegend = FALSE) %>%
  add_trace(x =~sq_lit, y =~LB, fill = 'tonexty', fillcolor='rgba(0,100,80,0.2)', 
            line = list(color = 'transparent'), name = 'Lower CI', showlegend = FALSE) %>%
  layout(xaxis = list(dtick = 10))

Interpretation:

• When the adult literacy rate is 67, a 1 point increase in inequality is associated with a 0.07 increase in life satisfaction, on average, controlling for population and GDP.

10.1 Zero Conditional Mean

A violation of the zero conditional mean assumption risks biased coefficients. We check for violations of the zero conditional mean in the plot below, which depicts our residuals by a linear combination of our independent variables. In the plot, we look for a roughly flat red line along zero. We can see that the red line is indeed mostly flat with some slight curvatures towards the ends, likely owing to the fact that we have limited observations. We therefore assess that our primary model does not violate this assumption.

plot(fit1, which = 1)

10.2 Homoskedasticity

The assumption of homoskedastic error consists of the following two components:

1. \(\epsilon_{i}\) is independent of \(\epsilon_{j}\) conditional on X for all \(i\) \(\neq\) \(j\).

2. The variance of the error does not depend on the predictor: Var(\(\epsilon_{i}\)|X) = Var(\(\epsilon_{i}\))

The assumption that all the errors have the same variance is called homoskedasticity.

The null for the Breusch-Pagan test is homoskedasticity. Given the results of our test; we reject the null and assume we have hetereoskedasticity. However, we need to be cautious about this test as it is sensitive to the number of observations.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(sandwich)
bptest(fit1)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit1
## BP = 17.795, df = 5, p-value = 0.003214

We might also suspect evidence of heteroskedasticity in our plot below which shows some unevenness in the scatter plot particularly around the center values as they cluster further toward the right. However, it does not look too bad.

plot(fit1, which = 1)

We can also judge heteroskedasticity from the next plot which if it showed a non-flat line. If we had homoskedastic data, the red line would be roughly flat, which is what we see.

plot(fit1, which = 3)

Taken collectively, we judge that we have homoskedasticity and this means that practically, we do not have to incorporate robust standard errors, although we still might want to in order to be conservative.

10.3 Normality of Errors

The normality of errors assumes that the errors must have this specific distribution:

\[\begin{eqnarray} \epsilon \sim N(0, \sigma^2), &\\ Corr(\epsilon_{j}, \epsilon_{k}) = 0, \forall j \neq k \end{eqnarray}\]

Which means:

1. all errors are normally distributed with a mean of 0,

2. all errors are independent from one another,

3. and all errors have the same variance.

If the errors are not normally distributed, or do not have a mean of 0, then we say we have a model misspecification problem. If the errors are not independent from one another, we have an autocorrelation problem.

We conduct several tests below which taken collectively, suggest that our errors are normal, as evidenced by: (1) the histogram of residuals which show a near normal distribution, (2) the shapiro test which fails to reject the null hypothesis that our residuals are drawn from a normal distribution, and (3) a normal Q-Q plot which shows only a few observations that deviate from the diagonal line at each end.

hist(fit1$residuals, breaks = 50, main = "Model 2 Residuals", xlab="Residuals")

shapiro.test(fit1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.98825, p-value = 0.5124

# qq plot
plot(fit1, which = 2)