Research Methods in Political Science I

Histograms

Let’s make a histogram of the variable height.

hist_height <- ggplot(myd, aes(height))
hist_height + geom_histogram()

## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

To save the figure created by ggpot(), we use ggsave() function. Fo instance,

ggsave('height-histogram.pdf')

saves the last figure you made as “height-histogram.pdf” in the working directory. In RStudio, you can save the figure in the bottom right panel by clicking “Export” on top fo the figure.

You can also draw graphs without ggplot2.

hist(myd$height, xlab = "height", ylab = "count", main = "", col = "gray")

What are diffreces between the two histograms we created?

The vertical axis of histograms is row count (or frequency) by default. Let’s use density for the vertical axis, instead.

hist_height + 
    geom_histogram(aes(y = ..density..))  

## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

Without ggplot2,

hist(myd$height, probability = TRUE, col = "gray",
     xlab = "height", ylab = "density", main = "")

The bin width (or the number of classes) changes the apperance of histograms, and we should manually choose the best bin width. Let’s change the bin width to 5.

hist_height + 
    geom_histogram(binwidth = 5)

Now, the histogram drawn with ggplot2 looks very similar to one without.

So far, we have treated female and male together. Let’s draw a histogram for each gender.

hist_height + 
    geom_histogram(binwidth = 5) + 
    facet_wrap( ~ sex)

As we can see, the range of height differs between gender. Let’s use an appropriate range for each sex.

hist_height + 
    geom_histogram(binwidth = 5) + 
    facet_wrap( ~ sex, scales = "free_x")

Box[-and Whisker] Plots

Next, let’s make box plots for female and male heights.

box_height <- ggplot(myd, aes(sex, height))
box_height + geom_boxplot() 

Scatter Plots

Now let’s examine the relationship between height and weight by a scatter plot.

scat_ht_wt <- ggplot(myd, aes(height, weight))
scat_ht_wt + geom_point()

Use different colors for different genders.

scat_ht_wt + geom_point(aes(color = sex))

Add a (regression) line to show a relationship between the variables.

last_plot() + geom_smooth(method = 'lm', se = FALSE, aes(color = sex))

Lastly, add income to the graph.

scat_ht_wt + geom_point(aes(color = sex, size = income))

Measures of Central Tendency

The mean of a variable \(x\) is denoted as \(\bar{x}\) (reads “x bar”): \[ \bar{x} = \frac{\sum_{i=1}^n x_i}{n}.\] In R, you can calculate the mean as follows.

x <- rnorm(20, mean = 10, sd = 4)  ## randomly draw 20 x's from N(10, 16)
sum(x) / length(x)

## [1] 10.28859

Since you will calculate the mean over and over for a lot of different variables, you might want to have a function to get the means. In R, you can create your own functions. Let’s make a function to obtain the mean.

get_mean <- function(x) {
    return(sum(x) / length(x))
}

Now you can get the mean by the function get_mean().

get_mean(x)

## [1] 10.28859

However, such a basic function is built in. You can get the mean by mean() (how easy!).

mean(x)

## [1] 10.28859

Similarly, you can obtain the median, the minimum, and the maximum by median(), min(), and max(), respectively.

median(x)

## [1] 10.64403

min(x)

## [1] 3.541893

max(x)

## [1] 14.40589

Quantiles

You can calculate quantiles by quantile(). You also use this function to get percentiles or quartiles.

## 50% quantile = 50th percentile = median
quantile(x, 0.5)

##      50% 
## 10.64403

## 25th percentile = 1st quartile and 75th pecentile = 3rd quartile 
quantile(x, c(.25, .75))

##       25%       75% 
##  8.739724 11.712385

Thus, you get the inter-quartile range (IQR) as follows.

quantile(x, .75) - quantile(x, .25)

##      75% 
## 2.972662

Alternatively, you can use

IQR(x)

## [1] 2.972662

The five number (0th to 4th quartiles) summary can be obtained by

quantile(x, c(0, .25, .5, .75, 1))   ## Five number summary with quartiles

##        0%       25%       50%       75%      100% 
##  3.541893  8.739724 10.644034 11.712385 14.405890

If you would like to use hinges instead of the 1st and 3rd quartiles (See Tukey 1977. Exploratory Data Analysis.), use fivenum().

fivenum(x)   ## Five number summary with hinges

## [1]  3.541893  8.711931 10.644034 11.820155 14.405890

Variance and Standard Deviation

The unbiased variance of a variance \(x\) is denoted as \(u_x^2\): \[u_x^2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n-1}.\]

In R, you can calculate the variance as follows.

sum((x - mean(x))^2) / (length(x) - 1)

## [1] 7.767437

The variance is the most important statistic, and the function to get it is of course built in. You can get it by var().

var(x)

## [1] 7.767437

The standard deviation is the square root of the variance:

sqrt(var(x))

## [1] 2.787012

You can also use the built-in function sd().

sd(x)

## [1] 2.787012

Population Data

If the vector \(x\) is not a sample but the population, the mean and the variance of x are denoted as \(\mu\) and \(\sigma^2\), respectively. The mean of the population is obtained by the same function as the sample mean.

mean(x)

## [1] 10.28859

However, the formula for the population variance is different than the sample variance: \[\sigma^2 = \frac{\sum_{i=1}^N (x_i - \mu)^2}{N}.\] Therefore, you cannot use var(). Instead, you should write

var(x) * (length(x) - 1) / (length(x))

## [1] 7.379065

There is no built-in function for this. Let’s make a function.

pop_var <- function(x) {
    return(sum((x - mean(x))^2) / length(x))
}
pop_var(x)

## [1] 7.379065

Simulation 1: Sample Means of Discrete Uniform

There is a bag containing ten balls numbered 0 through 9. If you randomly choose a ball from the bag, there are ten possible numbers you get, and they are equally likely. Thus, the distribution of the number you get by drawing a ball is discrete uniform, and the mean of the possible values is 4.5.

Now, suppose that you do not know the numbers written on the balls. You will infer the true mean (popoulation mean) of the numbers by the sample mean. You will obtain the sample mean by drawing a ball \(n\) times from the bag. You will draw only a ball at a time, and put it back to the bag before you draw the next ball. That is, you conduct sampling with replacemnet. We will run a simulation to examine the sampling distribution obtained by this procedure.

First, make a bag contaning ten balls numbered 0 through 9.

bag <- 0:9

Next, set the number of samples and the sample size (Note for Japanese speakers: 標本数と標本サイズの違いに注意).

trial <- 1000      ## the number of samples
n <- 1             ## sample size

Then, create a vector to save the simulation result. The length of the vector should match the number of samples. Before the simulation, we set each element to a missing value NA. (You can use any value to prepare a vector, but you should use NA because you can detect an error if the vector has NA after the simulation is run. If you use a numerical value, you cannot tell if it is the result of the simulation or the initial value.)

sample_means <- rep(NA, trial)

Now let’s run the simulation. We draw a ball from the bag with the function sample(). To conduct smapling with replacement, we set the argument relace = TRUE. Every time we sample \(n\) values, we calculate the sample mean, and save it to the vector sample_means. We repeat this for each sample using for loop. R code:

for (i in 1:trial) {
    balls <- sample(bag, size = n, replace = TRUE)
    sample_means[i] <- mean(balls)
}

Now sample_means contains the simulation result. Let’s draw a histogram to see a sampling distribution (Note that this is merely a pseudo sampling distribution since the true sampling distribution must have every possible sample and usually only a theoretical distribution).

hist(sample_means, prob = TRUE, col = "gray", axes = FALSE,
     breaks = 0:10, right = FALSE,
     xlab = "sample mean", ylab = "",
     main = "Sampling distribution from Uniform: n = 1")
axis(1, 0.5:9.5, labels = 0:9)
abline(v = 5, col = "red", lwd = 2)  ## mark the population mean

Does this look like normal? Probably not, because the sample size \(n\) is small (or not large enough). The histogram looks like the original uniform distribution, because \(n=1\).

What will happen as \(n\) gets larger? Run simulations! (See Assignment 2)