## 3.2 Solutions

### 3.2.1 Exercise 1

Create a new file called `assignment3.R`

in your `PUBL0055`

folder and write all the solutions in it.

In RStudio, go to the menu and select **File** > **New File** > **R Script**

Make sure to clear the environment and set the working directory.

```
rm(list = ls())
setwd("~/PUBL0055")
```

Go to the menu and select **File** > **Save** and name it `assignment3.R`

### 3.2.2 Exercise 2

Turn former colonies into a factor variable and choose appropriate labels.

First, we load the dataset.

`world.data <- read.csv("QoG2012.csv")`

Let’s take quick look at `former_col`

variable

`table(world.data$former_col)`

```
0 1
72 122
```

Use the `factor()`

function to convert `former_col`

to a factor variable. According to the codebook, 1 means colonised and 0 means not colonised.

```
world.data$former_col <- factor(
world.data$former_col,
labels = c("not colonised", "colonised"),
levels = c(0,1)
)
```

### 3.2.3 Exercise 3

How many countries were former colonies? How many were not?

We can get the numbers from a frequency table.

`table(world.data$former_col)`

```
not colonised colonised
72 122
```

122 countries were colonised

You can also count them by each category individually

`sum(world.data$former_col == "colonised")`

`[1] 122`

`sum(world.data$former_col == "not colonised")`

`[1] 72`

### 3.2.4 Exercise 4

Find the means of political stability in countries that (1) were former colonies, (2) were not former colonies.

We use the `mean()`

function to get the mean of political stability and subset for the groups of countries using the square brackets.

`mean(world.data$wbgi_pse[world.data$former_col=="colonised"])`

`[1] -0.231612`

`mean(world.data$wbgi_pse[world.data$former_col=="not colonised"])`

`[1] 0.2858409`

The average level of political stability in countries that were not colonised is 0.29. Mean political stability in countries that were colonised is -0.23.

The variable political stability `wbgi_pse`

is an index. Larger values correspond with more political stability. We see that political stability is higher in countries that were not colonised.

Looking at this difference, we might conclude that the legacy of colonialism is still visible today and manifests itself in lower political stability. Let’s investigate further to see whether the difference in means is statistically significant.

### 3.2.5 Exercise 5

Is the difference in means statistically significant?

We can test whether the difference is statistically significant using the t-test.

`t.test(world.data$wbgi_pse ~ world.data$former_col, mu = 0, alt = "two.sided")`

```
Welch Two Sample t-test
data: world.data$wbgi_pse by world.data$former_col
t = 3.4674, df = 139.35, p-value = 0.0006992
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.2224004 0.8125053
sample estimates:
mean in group not colonised mean in group colonised
0.2858409 -0.2316120
```

The p-value is small. Smaller than the conventional alpha level of 0.05. We can also look at the confidence interval which ranges from 0.22 to 0.81. So, if we were to repeatedly sample, the confidence interval of each sample would include the true population mean \(95\%\) of the time. Or more intuitively, we are \(95\%\) confident that the average population level of political stability is within our interval.

Let’s look at the variable a little closer to understand this difference in substantive terms.

First we compute the difference in means.

```
mean(world.data$wbgi_pse[world.data$former_col=="not colonised"]) -
mean(world.data$wbgi_pse[world.data$former_col=="colonised"])
```

`[1] 0.5174529`

The numerical difference between the two means is 0.52. Is this difference small or large? That is difficult to say because the variable is not measured in intuitive units (income in dollars is an example of a variable that is measured in intuitive units).

Let’s look at the variable a little closer to understand this difference in substantive terms.

`min(world.data$wbgi_pse)`

`[1] -2.467461`

`max(world.data$wbgi_pse)`

`[1] 1.675609`

So the range of the variable is from -2.47 to 1.68

What does the difference of 0.52 mean in terms of the range of values for this variable? We can divide this difference by the range of the variable to better understand what it means in relative terms:

\[ \frac{0.517}{1.676 - (-2.467)} \]

Getting the difference as a percentge of the range is helpful to understand the variable in substantive terms. The difference in means is 0.52. That constitutes 0.12 (12.49%) of the range of the variable. So the difference in political stability between countries that were not colonised and those that were colonised is indeed large.

### 3.2.6 Exercise 6

In layman’s terms, are countries which were former colonies more or less stable than those that were not?

The results from the t-test show that countries that were colonies are less stable than those that were not. The difference is large and statistically significant.

### 3.2.7 Exercise 7

How about if we choose an alpha level of 0.01?

The p-value is 0.000699. That is smaller than an alpha level of 0.01 as well. Therefore, picking an alpha level of 0.01 does not change our results.

### 3.2.8 Exercise 8

What is the level of measurement of the United Nations Human Development Index (HDI) variable `undp_hdi`

?

Acccording to the United Nations Development Programme, Human Development Index(HDI) is a composite index of three key dimensions: life expectancy, education, and per capita income. The description goes on to tell us how the dimensions are measured and that the index is the geometric average of the components. The variable is, therefore, continuous.

### 3.2.9 Exercise 9

Check the claim that its true population mean is 0.85.

Let’s estimate the mean from our sample

```
hdi_mean <- mean(world.data$undp_hdi, na.rm = TRUE)
hdi_mean
```

`[1] 0.69824`

Our estimate is 0.70. The claim is that it is 0.85.

**Null hypothesis**: The true population mean of the human development index is: 0.85.**Alternative hypothesis**: The true population mean is different from 0.85.

We pick an alpha level of 0.05 for our test and use 0.85 as the value for the `mu`

argument.

`t.test(world.data$undp_hdi, mu = 0.85, alt = "two.sided")`

```
One Sample t-test
data: world.data$undp_hdi
t = -11.139, df = 174, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0.85
95 percent confidence interval:
0.6713502 0.7251298
sample estimates:
mean of x
0.69824
```

`t_test_result <- t.test(world.data$undp_hdi, mu = 0.85, alt = "two.sided")`

The p-value is lower than 0.05 and hence we reject the null hypothesis (that hdi is 0.85). Looking at our confidence interval, we expect that if we were to repeatedly sample, the population mean would fall into the interval 0.67 to 0.73, \(95\%\) of the time.

**Scientific Notation**

The p-value above includes an \(e\) symbol that some of you may not have seen before. The \(e\) symbol is used in scientific notation to represent very large or very small numbers.

The \(e\) in the p-value 2.2e-16 stands for \(\times 10^{-16}\). So the p-value is:

\[ \begin{split} p & = 2.2e^{-16} \\ & = 2.2 \times 10^{-16} \\ & = \frac{2.2}{10^{16}} \\ & = \frac{2.2}{10000000000000000} \\ & = 0.00000000000000022 \end{split} \]

There’s an excellent video tutorial from Khan Academy if you need additional help with scientific notations.

### 3.2.10 Exercise 10

Calculate the t-statistic.

We could take the t-statistic from the t-test above. But we are going through the individual steps here.

Step 1: We get the number of non-missing observations of `undp_hdi`

.

`n <- sum(!is.na(world.data$undp_hdi))`

Step 2: We calculate the standard variation of `undp_hdi`

.

`hdi_sd <- sd(world.data$undp_hdi, na.rm = TRUE)`

Step 3: We calculate the standard error of the mean which is:

\[ \sigma_{\bar{Y}} = \frac{ \sigma_{Y}} { \sqrt{n}} \]

`se.y_bar <- hdi_sd / sqrt(n)`

Step 4: We compute the t-statistic as the difference between our estimated mean and the mean under null and then divide by the standard error of the estimated mean.

\[ t = \frac{ \bar{Y} - \mu_{0}} { \sigma_{\bar{Y}} } \]

Note: By dividing by the standard error of the estimated mean, we are making the difference of the means (our estimated mean - the mean under the null hypothesis), free of the units that the variables were measured in. The t-statistic is the difference of the means, where its units are standard deviations.

The t-statistic follows a t-distribution with n-1 degrees of freedom. It is well approximated by a standard normal distribution.

```
t.statistic <- (hdi_mean - 0.85) / se.y_bar
t.statistic
```

`[1] -11.13908`

The t.statistic is -11.14. Given that this t-statistic follows a standard normal distribution approximately (because the sample is large enough), we would reject the null hypothesis if the t-statistic is more than 1.96 standard deviations away from 0. This t-statistic is more than 11 standard deviations away from zero. The t-statistic is very far in the tails of the distribution. It is very unlikely that we would have estimated the mean that we did (0.70) if the population average of HDI really was 0.85.

Therefore, we can reject the null hypothesis.

### 3.2.11 Exercise 11

Calculate the p-value.

We have the t-statistic already. We will follow the steps to get the p-value. But first, take a look at what the distribution of our t-statistic looks like. Our t-statistic is -11.14. This is very very far in the left tail.

Step 1: We compute the cumulative probability to get our t-statistic or something smaller from this distribution.

```
cumulative.prob <- pt(t.statistic, df = n-1)
cumulative.prob
```

`[1] 2.127538e-22`

Notice that our t-statistic is negative. By calculating the probability to get a t-statistic of -11.14 or smaller, we get the proability that we would see a t-statistic like this in the left tail.

Step 2: Since we are conducting a two-sided test, we want the proability in the right tail as well.

In other words, we want the probability of getting a t-statistic that is 11.14 or larger. The distribution is symmetric so in the next step we can just multiply the probability by 2.

`cumulative.prob * 2`

`[1] 4.255075e-22`

The p-value is extremely small. In fact, it is so small that the t-test did not report the correct p-value either but simply stated that the p-value is smaller than 2.2e-16.

Notice also that in the seminar we calculated `(1 - pt(t.statistic, df = n-1)) * 2`

. This was because our t-statistic was positive in the seminar. Here it is negative so we don’t need to subtract the cumulative probability from 1.

### 3.2.12 Exercise 12

Construct a confidence interval around your mean estimate.

There are two ways to do this which yield slighly different answers.

**Approach 1**

Our sample is large. The t-statistic follows approximately a standard normal distribution. We can, therefore, use the critical value from the standard normal distribution. \(95\%\) of the area under the standard normal distribution is within \(1.96\) standard deviations from the mean.

The lower bound for the confidence interval is:

`hdi_mean - (1.96 * se.y_bar)`

`[1] 0.6715367`

Similarly, the upper bound is:

`hdi_mean + (1.96 * se.y_bar)`

`[1] 0.7249432`

Our best guess of the population mean is 0.70. We are \(95\%\) confident that our interval from 0.67 to 0.72 includes the population mean.

**Approach 2**

We get the critical value from a t-distribution with n-1 degrees of freedom. To find the critical value, we use the t-distribution’s quantile function `qt()`

.

```
t.critical <- qt(0.975, df = n-1)
t.critical
```

`[1] 1.973691`

So, in our t-distribution \(95\%\) are within precisely 1.97369143975607 standard deviations from the mean.

Notice, that the confidence interval covers \(95\%\) of the area under the curve and is centered on the mean. There are \(2.5\%\) in each tail. `qt()`

gives us the critical value at the cumulative probability of whatever we enter into the function. We input 0.975, so we get the critical value at the point where the right tail starts.

`hdi_mean - (t.critical * se.y_bar)`

`[1] 0.6713502`

`hdi_mean + (t.critical * se.y_bar)`

`[1] 0.7251298`

Using this more exact approach, we show that the \(95\%\) confidence interval covers values of the human development index from 0.67 to 0.73.

### 3.2.13 Exercise 13

Discuss your findings in terms of the original claim. Interpret the t-statistic, the p-value, and the confidence interval.

We can reject the original claim that the mean of HDI is 0.85. Given, an alpha level of 0.05, our large t-statistic (-11.14) implies that our estimate (0.70) is extremely unlikely assuming that 0.85 is the population mean. From exercise 11, we know the p-value which is the probability that we have mistakenly rejected a correct null hypothesis. That probability is 0.00

Our best guess of the population mean is between 0.67 and 0.73 (rounded to two digits).

### 3.2.14 Optional Exercises that require reading Extra Info below

#### 3.2.14.1 Optional Exercise 1

Create a scatter plot with latitute on the x-axis and political stability on the y-axis.

The relationship between two continuous variables is best illustrated with a scatterplot.

```
plot(
wbgi_pse ~ lp_lat_abst,
data = world.data,
pch = 19,
frame.plot = FALSE,
xlim = c(0, 0.9),
main = "Relation between pol. stability and distance from the equator",
xlab = "Distance from equator in decimal degrees (0 to .9)",
ylab = "Political stability index"
)
```

#### 3.2.14.2 Optional Exercise 2

What is the correlation coefficient of political stability and latitude?

`cor(world.data$wbgi_pse, world.data$lp_lat_abst, use = "complete.obs")`

`[1] 0.3936597`

the correlation coefficient ranges from -1 to 1 and is a measure of linear association of two continuous variables. The variables are positively related. That means, as we move away from the equator, we tend to observe higher levels of political stability.

#### 3.2.14.3 Optional Exercise 3

If we move away from the equator, how does political stability change?

Political stability tends to increase as we move away from the equator.

#### 3.2.14.4 Optional Exercise 4

Does it matter whether we go north or south from the equator?

It does not matter whether we go north or south. Our latitude is measured in decimal degrees. Thus, a value of 0.9 could correspond to either the North Pole or the South Pole.

### 3.2.15 Advanced Exercises

#### 3.2.15.1 Advanced Exercise 1

Calculate the numerical difference in means (political stability conditional on colonialization) using the `mean()`

function.

We actually solved this as part of Exerise 5 above.

```
mean(world.data$wbgi_pse[world.data$former_col=="not colonised"]) -
mean(world.data$wbgi_pse[world.data$former_col=="colonised"])
```

`[1] 0.5174529`

The difference in means is 0.52. Countries that were not colonised are more politically stable.

#### 3.2.15.2 Advanced Exercise 2

Calculate the standard error of the difference in means (hint: using just the `sd()`

function is incorrect in this context).

The standard error of the difference between two means is:

\[ \sqrt{\frac{\sigma^2_{Y_A}}{n_A} + \frac{\sigma^2_{Y_B}}{n_B}} \]

Where \(A\) and \(B\) are the two groups. In our case the two groups are colonised and not colonised:

\[ \sqrt{\frac{\sigma^2_{Y_{colonised}}}{n_{colonised}} + \frac{\sigma^2_{Y_{not.colonised}}}{n_{not.colonised}}} \]

```
colonised <- world.data[world.data$former_col=="colonised", ]
not.colonised <- world.data[world.data$former_col=="not colonised", ]
sqrt(
(var(colonised$wbgi_pse) / nrow(colonised)) +
(var(not.colonised$wbgi_pse) / nrow(not.colonised))
)
```

`[1] 0.1492324`

The standard error of the difference in means is 0.15.

#### 3.2.15.3 Advanced Exercise 3

Is the difference in means more than \(1.96\) standard deviations away from zero? Interpret the result.

0.5174529 - (2 * 0.1492324) = 0.22

The difference in means is further than two standard deviations away from zero. Given an alpha level of 0.05, we can reject the null hypothesis that political stability is the same in countries that were colonised and those that were not.

#### 3.2.15.4 Advanced Exercise 4

We claim, the difference in means in terms of political stability between countries that were former colonies and those that were not is 0.3. Check this hypothesis.

We use the t-test for the difference in means to test this claim. Normally, our null would be that the difference in means is zero, i.e. there is no difference. Here, the claim is that the difference is 0.3. So, all we have to do is adjust the null hypothesis accordingly.

`t.test(world.data$wbgi_pse ~ world.data$former_col, mu = 0.3, alt = "two.sided")`

```
Welch Two Sample t-test
data: world.data$wbgi_pse by world.data$former_col
t = 1.4571, df = 139.35, p-value = 0.1473
alternative hypothesis: true difference in means is not equal to 0.3
95 percent confidence interval:
0.2224004 0.8125053
sample estimates:
mean in group not colonised mean in group colonised
0.2858409 -0.2316120
```

We cannot reject the claim that the true difference in means is 0.3.

#### 3.2.15.5 Advanced Exercise 5

An angry citizen who wants to defund the Department of International Development (DFID) claims that countries that were former colonies have reached \(75\%\) of the level of wealth of countries that were not colonised. Check this claim.

The null hypothesis is that there is no difference between the level of wealth in countries that were former colonies and 0.75 times the level of wealth in countries that were not former colonies.

The claim of the angry citizen is:

`mean(not.colonised$wdi_gdpc, na.rm = TRUE) * 0.75 `

`[1] 12311.54`

Estimated level of wealth in countries that were colonised

`mean(colonised$wdi_gdpc, na.rm = TRUE)`

`[1] 6599.714`

Our estimate is roughly half the angry citizen’s claim. The substantial difference is huge. To cover all our bases, let’s perform the t-test for the difference in means.

How would we do this though? The trick here is to provide vectors of GDP values for each group (colonised and not-colonised) as two arguments to the `t.test()`

function instead of using a formula.

`t.test(not.colonised$wdi_gdpc * 0.75, colonised$wdi_gdpc, mu = 0, alt = "two.sided")`

```
Welch Two Sample t-test
data: not.colonised$wdi_gdpc * 0.75 and colonised$wdi_gdpc
t = 3.6218, df = 127.72, p-value = 0.0004206
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
2591.29 8832.37
sample estimates:
mean of x mean of y
12311.544 6599.714
```

Clearly, we can reject the citizen’s claim. Our p-value implies that the probability that we see this huge difference in our data, given that there really is no difference, is much smaller than our conventional alpha level of 0.05.